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1. Preface 

That ^^differentiating is an operation easier than integrating^^ is a statement that should not sound 
too surprising; while, a more pleasant wonder might result when suitable differentiations make us 
reduce, if not avoid at all, the number of direct integrations - of course the two operations, being 
the inverse of each other, have not to be thought as performed with respect to the same variable! 
As paradigmatic example, let us just consider the class of integrals, 

In{a) = [ e""^' a;" dx . 
Jo 

For n = 0, this is just the Gaussian integral, 

while for n = 1, the integrand is integrable by quadrature. 



h{a) = 




To compute /„ with n > 1, one can use the identity, 

^ I - 
oa 

In fact, 

• n = 2s 




• n = 2s + l 

Therefore the infinite set of integrals /„ can be computed without any integration, provided the 
knowledge of just two basic integrals, namely Jo, and Ii, that in the forthcoming terminology would 
be defined as the master integrals of the class /„. 

The above example was a too lucky one: i) the repeated a-derivative did not entangle integrals 
having even and odd indices, therefore Jg and Ji never appear linked by any differential identity; ii) 
the value of the master integrals was known, possibly obtained by direct integrations. 

In the more general case masters' are unknown, and their evaluation becomes an open problem. In 
the following pages, we will sec how the exploitation of integration-by-parts not only yields algebraic 
relations among infinite sets of integrals and their masters', but as well leads to differential equations 
satisfied by the master integrals themselves. Solving these differential equations becomes a tool for 
computing master integrals, when their direct integration is not viable. 
As it happens to (many) Feynman integrals. 
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2. Introduction 

A perturbative approach to the quantitative description of the scattering of particles in quantum 
field theory involves the computation of Feynman diagrams. For a given number of external particles 
- the legs of diagram - fixed by the process under study, and a given order in perturbation theory, 
the skeletons of diagrams are built up by joining the edges of legs and propagators into vertexes, 
forming tree patterns and closed loops. 

Beyond the tree level, each Feynman diagram represents an integral which has, in general, a tenso- 
rial structure, induced by the tensorial nature of the interacting fields. Therefore, the result of its 
evaluation must be a linear combination of the tensors provided by the theory and by the kinematics 
of the process under study. The coefhcients of this linear combination, usually called form factors, 
can be always extracted from each Feynman diagram, before performing any evaluation, by means 
of suitably chosen projectors. 

These form factors are scalar integrals closely connected to the original Feynman diagram: the nu- 
merator of their integrand may contain all the possible scalar products formed by external momenta 
and loop variables; whereas its denominator is formed by the denominators of propagators present 
in the diagram itself. 

Due to the bad convergence of loop integrals in four dimensions, regularization prescriptions 
are mandatory. Hereafter the integrals are r egul arised within the framework of 't Hooft-Veltman 
continuous-dimensional regularisation schemJi2Sl_ Accordingly, the dimension D of an extended in- 
tegration space is used as a regulator for both infrared ( IR) and ultraviolet (UV) divergences, which 
finally do appear as poles in {D — 4) when D goes to 4 USSI 

The aim of a precise calculation is to compute Feynman diagrams for any value of the available kine- 
matic invariants. Except in case of simple configurations (e.g. very few legs and/or few scales), quite 
generally, approximations have to be taken by limiting the result to a specific kinematics domain, 
and, thus, looking for a hierarchy among the scales, to get rid of the ones which anyhow would give 
a negligible contribution in that domain. 

The puzzling complexity of the Feynman diagrams calculation arises because of two different 
sources: either multi-leg or multi-loop processes. In recent years the progress in the evaluation of 
higher loop radiative corrections in quantum field theory has re ceived a strong boost, due to the 
optimisation and automatising of various techniques (see refs. in I 59 | 79 [77|| ^-j^jg -^^Qp]^ -^g review 
one of the most effective computational tools which have been developed in the framework of the 
dimensional regularization: the method of differential equations for Feynman integrals. 

The method was first proposed by Kotiko^ffi in the early nineties, while dealing with the evaluation 
of 2- and 3-point functions. The basic idea was to consider a given unknown integral as a function 
of one of the propagator masses, and to write for it a differential equation in that variable. Thus, 
instead of addressing its direct integration, the value of the integral could be found by solving the 
differential equation. 

The advantages of the novel ideas were soon realized I^EEEEI^ g^j^^j generalised at a later stage by 
RemiddilS, who proposed the differentiation with respect to any other available kinematics invariants 
formed by the external momenta. That enabled the application of the differential equation method 
also to integral with massless propagators (provided the existence of any other non-trivial scale). 
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Finally, Remiddi and Gehrmann l ^ l ^ l ^ l fully developed the method by showing its effectiveness 
through the systematic application to a non-trivial class of two-loop four-point functions, whose 
result is still considered as state-of-the-art. 

From that mome nt on, ,j'.^6,_™6thod of differentia l equations became to be widely used 
in different contexts | 11 | 12 | 13 | 14|15|16 | 17 | 18|19 | 20 | 21|22 | -pj^^ jj^^^ unprecedented results ob- 
tained thr ough its application spans among multi-loop functions from zero to four ex- 
ternal legs, I^5l27l28l^9l5^6^^3ll3:jl33l3ll35l36l37l38r39l40l41l4:j|43l44l4^ 

and within the most interesting sectors of particle phenomenology, like jets physics 
[8 9 10 50 51 52 53 54 55 M 5Z 58 correct ions to lepton form factors ' ^"'^^l^^ i Bhabha 

Scattering I 109|11 0|l ll | 112 | 113 | 114 i ll5 | 116 | 117 |. q^d corrections to lept on form facto rs and 
top-phy sics t^ '^^^ 85 86 87 1^ forw ard-backward asym metry of heavy-quark I 9 5 | 96 | 97 | 98 | 99l . Higgs 
Physics 1 05 106.10^ , Electrowcak sector | 121 | 122 | 123 | 124 | 125 | 126 [ g^^akov form fac- 

tors 1^^^ se mileptonic decay l^l^ ^l^-l^ static parameters and gauge boson properties 

l88l89l90l9ll92ig5iMl 

The efforts to achieve analytic solutions of differential equations for Feynman integrals has stimu- 

1 97 1 98 

lated new developments on the more mathemati cal side .^^"i^ especially concerning the properties 
of transcendental functionsl l^'^ l l^^ l l^^ l l^^ l l^S I jj^ particular, a novel s et of function s that generalize 
Nielsen polylogarithms, the so called Harmonic Polylogarithms (HPL) | 129 | 130 |131J^ have been found 
suitable for casting the result in analytic form - that means without ambiguities due to zeroes hid- 
den in functional relations, and supplied with series expansions. While HPL's can be considered as 
iterative integral of rational kernels, recently it has been pointed out that the solution of differential 
equations for generic inte grals wit h massive loops demands as well for irrational integration kernels, 
yielding elliptic functionslI3HI102] 

The range of applicability of the differential equations technique is broadened by the possibility 
of a natural switch toward a semi-numerical approach, since, whenever the analytic integration were 
not required or not viable, the differential equation(s), analytically obtained, could be solved with 
numerical techniques . 

Nowadays we are not at the point t o have the method for eyaluating any Feymnan integral, 
but certainly we dispose of several tools I «" l «l | 71 l «^ l «5 | 63 | 64 | 67 | 68 | 72 | 73 | 69l70 | 74 | 75 | 76 | 77 | ^ attack 
successfully many problems in perturbation theory, and usually a combination of them is necessary 
for the achievement of the final answer. Therefore let us discuss in detail how to build ans solve 
differential equations for integral associated to Feynman graphs. 

The computational strategy is twofold. 

• In a preliminary stage, by exploiting some remarkable properties of the dimensionally regu- 
larised integrals, namely integration-by-parts identities (IBP), Lorentz invariance identities 
(LI), and further sets of identities due to kinematic symmetry specific of each diagram, one 
establishes several relations among the whole set of scalar integrals associated to the original 
Feynman diagram. 

By doing so, one reduces the result, initially demanding for a large number of scalar inte- 
grals (from hundreds to billions, according to the case), to a combination of a minimal set 
(usually of the order of tens) of independent functions, the so called master integrals (Mi's). 



Feynman Diagrams & Differential Equations 5 



• The second phase consists of the actual evaluation of the Mi's. By using the set of identities 
previously obtained, it is also possible to write Differential Equations in the kinematic 
invariants which are satisfied by the Mi's themselves. When possible, these equations can 
be solved exactly in D dimensions. Alternatively, they can be Laurent-expanded around 
suitable values of the dimensional parameter up to the required order, obtaining a system 
of chained differential equations for the coefhcients of the expansions, which, in the most 
general case, are finally integrated by Euler's variation of constants method. 

One of the key advantages of the method is that it yields a clear separation between the merely 
algebraic part of the work - which is, not surprisingly, always very he avy i n multiloop calculations, 
and can be most conveniently processed by a computer algebra progranll^^i^ from the actual analytic 
issues of the problem, which can then be better investigated without the disturbance of the algebraic 
complexity. 

The paper is organised as follows. In section 2, it is described how to reduce a generic (combination 
of) Feynman integrals to a limited set of Mi's and how to write the system of differential equations 
they fulfil. As illustrative applications, respectively in section 3-4, the one-loop and two-loop vacuum 
polarization functions in QED are explicitly computed. In the further two sections, we discuss some 
less obvious example of differential equations. In section 5, we describe the solution of a system of 
three coupled first-order differential equations, to compute three Mi's associated to a class of two- 
loop 3-point functions, addressing as well the problem of finding their boundary conditions. While 
in section 5, we describe the solution of a fourth-order differential equation involving the Mi's of a 
4-loop 2-point diagram, and it will be considered the link between differential and difference equation 
for Feynman integrals'^. 

3. Reduction to Master Integrals 
3.1. Topologies and Integrals 

Let us consider a Feynman diagram with i loops, g external legs, d internal lines and a given 
tensorial structure. Such a diagram, when not representing a scalar quantity, can be decomposed as 
a combination of products of a scalar form factor times a tensor. Thus, computing the contribution 
of any diagram to a given process is equivalent to the computation of scalar factors, which, after 
some preliminary algebra (evaluation of Dirac traces and contraction of Lorentz indices) reads as a 
combination of scalar integrals like. 



• Si (i — 1, • • • , Ngp) is any scalar product formed by either one of independent external 
momentum and an internal loop momentum, or by two internal loop momenta [ui is an 
integer exponent such that ni > 0) and N^p is the total number of such scalar products. 




d^h nsgf' 



(1) 



where 



6 M. Argeri & P. Mastrolia 



given by 

s.p. external — internal " 

s.p. internal — internal 

• "Dj — Qj + ni?j {j — 1, ■ • • ,d) is the denominator of the j-th propagator being qj and rrij, 
respectively, the corresponding momentum and mass. From now, we call Vj propagator. 

The first task relies on a suitable classification of the integrals, in order to minimize the number 
of the ones which have to be actually integrated. 

As a preliminary set up, one is invited to consider not anymore diagrams but topologies, that, 
drawn exactly as scalar Feynman diagrams, contain only and all different propagators (and scalar 
vertices). In view of settling down a correspondence between a given topology and the class of 
integrals it represents, one proceeds by simplifying as well the number of scalar products in the 
numerator. The number N^p is a redundant quantity and can be easily reduced, with a procedure 
commonly called trivial tensor reduction, according to which some (when not all) of the scalar 
products can be expressed in terms of the denominators of the topology. 

In general, if t is the number of propagators of a given topology, we can express t of the N^p 
scalar products containing the loop momenta in terms of the denominators - which will be later 
on simplified against the corresponding term in the denominator. Thus the most general integral 
associated to every topology does contain only the remaining q (= Ngp — t) irreducible scalar products 
and reads as 

d^fci d°k2 d^ke nLi^r' 

\D-2 (Oir\D-2 ' ' ' (0^\D-2 m -r^™! ' V'^-' 



(27r)^-2 (27r)^-2 (27r)^-2 Ii]^{D\^ 



where > and rrij > 1. 



One can denote with It^r.s the class of the integrals with: a given set of t denominators; q — N^p — t 
irreducible scalar products; a total of r = '^^ {mi — 1) powers of the t denominators; and s = nj 
powers of the q scalar products. It can be shown that the number of the integrals belonging to the 
class It,r,s is 

In the next sections, we will see that integrals of the type ([3]), belonging to a given topology, 
therefore differing for the values of the indices m^, n^, are not independent. Algebraic relations among 
them, can be written in the form of a sum of a finite number of terms set equal to zero, where each 
term is given by the product of a polynomial (of finite order and with integer coefficients in the 
variable D, masses and Mandelstam invariants) and one of the integrals belonging to It,r,s- They 
can be used recursively to express as many as possible integrals of a given class in terms of as 
few as possible (suitably chosen) ones. The way to generate those kind of relations goes through 
integration- by-parts, Lorentz invariance and symmetry considerations. 

Before going on with the discussion let us see explicitly how the general definition we have 
introduced do apply in practice. 
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3.1.1. Example: from the diagram to the topology 

Let us consider the Feynman diagram showed in Fig[T] 




ki — k 

Fig. 1. Feynman diagram with g=2 legs and 1 = 2 loops 

with 3 = 2 external legs, 1 = 2 loops and d = 5 internal lines, which gives a number of scalar products amounting to, 

Af^p = 2 (^2 + 1 - i 

The most general set of scalar integrals possibly arising from its computation has the following representation, 



(27r)0-2 (27r)D-2 Ii]^^Vj 

where 



(5) 





= g- 


h vn? 








Vz 


= (fci 


-k2) 


Da 


= (fci 






= X>i 


= kl- 



- m," 
^2 



5i 


= 1 




S2 


— 




53 


= fci 




54 


= k ■ 


ki 


S5 


= k ■ 


k2 



The original diagram contains five internal lines, but two propagators are indeed equal, so there are only four different 
propagators. Therefore, the integrals JSj actually belong to the simpler set. 



/ 



(27r)0-2 (27r)0-2 DlDiV^D^ 

which can be represented by the topology in Fig. |2] 

The trivial tensor reduction of the scalar product can be realized according to the following table, where one 
chooses to express four (out of five) scalar products in terms of the denominators. In the end, only one of the five 



Scalar product Si 


Corresponding propagator 


Relationship 


Si = kl 


Vi^kf + rri^ 


kl^Vi- rri^ 


S2 = kl 


V2 = kl 


kl^V2 


S3 = ki ■ k2 


V3 = {kl - k2f + 


kl-k2 = \{Vi^V2-V3) 


S4 = k ■ kl 


Vi = (fcl - + TO^ 


k-ki^ liVi-Vi + P) 



scalar products involving the loop momenta is left over as irreducible, arbitrary chosen to be fc2 • k. Therefore, the 
integrals in ((Sjl, represented by the topology in Fig|2] indeed belongs to the class, 

f d^ki d^ki (fc2-fc)"l 

The trivial tensor reduction might as well lead to the complete cancellation of some denominator. Should it be the 
case, the resulting integral can be classified as belonging to the subtopology obtained by pinching the corresponding 
internal line. 
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ki — k 

Fig. 2. 4-denominators topology 

3.2. Integration-by-parts Identities 

Integration-by-parts identities (IBP-Id's) are among the most remarkable properties of dimensionally 
regularized integrals and they were first proposed in the eighties by Chetyrkin and 

Tkachovll^a. The 

basic idea underlying IBP-Id's is an extension to D-dimensional spaces of Gauss' theorem. For each 
of the integrals defined in equation ([3]) one can write the vanishing of the integral of a divergence 
given by, 



In the above identities the index i runs over the number of loops {i — 1,2,...,£), and the vector 
can be any of the {£ + g — 1) independent vectors of the problem:fci, • • • , ki,pi, ■ ■ ■ ,Pg-i\ in such 
way, for each integrand, l{l + g — \) IBP-Id's can be established. When evaluating explicitly the 
derivatives, one obtains a combination of integrands with a total power of the irreducible scalar 
products equal to (s— 1), s and {s + 1) and total powers of the propagators in the denominator equal 
to {t + r) and (t + r + 1), therefore involving, besides the integrals of the class Ir,s,t, also the classes 
It,r,s-l^ It,r+i,s and It.r+i,s+i- Simplifications between reducible scalar products and propagators 
in the denominator may occur, lowering the powers of the propagators. During that simplification, 
some propagator might disappear, generating an integral belonging to a subtopology, with t — 1 
propagators. 



3.3. Lorentz Invariance Identities 

Another class of identities can be derived by exploiting a general properties of the integrals in ([4]), 
namely their nature as Lorentz scalars. If we consider an infinitesimal Lorentz transformation on 
the external momenta, Pi Pi + Spi, where Spi = iiJfivPi,v with a;^^ a totally antisymmetric tensor, 
we have 



I{p^+5pi)^I{p,). (9) 

Because of the antisymmetry of and because 

I[p, + 5p,) = I{p.) -f - ^(ft) + ^,uY.Pn/-P^^ (10) 

we can write the following relation 
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dp„ 



I{P^) = 0. 



(11) 



Eq. pip can be contracted with all possible antisymmetric combinations of the external momenta 
Pi,fiPj,ui to obtain other identities for the considered integrals. 

In case of integral associated to any vertex topologies with two independent external momenta, pi 
and p2 , we can build up the identity 

pi 



{Pi ■ P2) Pi 



d 



d 



d 

P2Pl.fJ.- 



d 



-P3 = 



(12) 



P2 



dpi,p dp2,p Op2.^i OPl,l-i, 

In the case of a richer kinematics, like in the case of integrals associated to box topologies with three 
independent external momenta, pi, p2 and ps, we ca write down three Ll-id's 



(pi,mP2,p - P1..uP2.j,)Y^ {pn,u 



d 



dpn„ 



d 



dpn,. 



Pi 



P2 



P3 



Pi 



(13) 



iPl,p.P3,f^-Pl,uP3,p.)^{pn,^- 



dpn„ 



Pn,ii 



dpn,v 



Pi 



P2 



P3 



P4 



(14) 



(P2,AiP3,A* - P2,uP3,p.) E ( 



' dp„ 



dpn,^ 



Pi 



P2 



P3 



= . 



P4 



(15) 



3.4. Symmetry relations 

In general, further identities among Feynman integrals can arise when it is possible to redefine the 
loop momenta so that the value of the integral itself does not change, but the integrand transforms 
into a combination of different integrands. By imposing the identity of the original integral to the 
combination of integrals resulting from the change of loop momenta, one obtains a set of identities 
relating integrals belonging to the same topology. 

More identities, originally found by Larin, may arise as well when there is a sub-loop diagram 
depending - after its integration - on a specific combination of momenta. By equating the original 
integral and t he o ne obtained by projecting onto such a combination of momenta, one gets additional 
relations (see for details). 



3.4.1. Example 1: IBP-id's 

Let us take again the two- loop topology of Fig. 2, and its class of integrals The corresponding IBP-Id's can be 
established by the vanishing of the following integral of a divergence, 



/ 



d^fci d°k2 d ( (fc • fc2)"i 



(i = l,2), (16) 



10 M. Argeri & P. Mastrolia 



where = ki,k2,k, which amount to 3 X 2 = 6 identities. 

For simplicity let us choose ni = and mi = ... = 1714 = 1. We have 

The six IBP-Id's are resumed in the following table 



derivative ttt— 


vector 


corresponding IBP-Id 


a 








= 








9 




f d^fei d°*:2 d 


( '=2,Ai 


= 




J (27t)"--^ (27r)"-^ 9fci,,, 


[^X'lX>2'D3l'4 J 


9 




r d^ki d°k2 d 




= 




J (27r)0-2 {2it)0-2 Qki_^ 


i •Dl-D2l'3f4 J 


d 




f d^ki d"k2 d 


ki,fi 


= 




J (27r)0-2 (27r)0-2 dk2,^ 


{ ViV2Vs,Vi ) 


a 




C d^ki d°k2 d 


( '=2,Ai 


= 


9fc2,f, 


J (27r)0-2 (27r)0-2 9fe2,M 


I X'lX)2X)3X)4 J 


d 








= 


9^2, f. 


J (27r)0-2 (2;r)0-2 9fe2,,, 


[■DlI52l'3X)4 J 



Fig. 3. IBP for the topology of FigO 



Let us perform explicitly the calculation of the last of above identities. 
We have 



/ 



d^ki d^k2 



a 







(27r)D-2 (27r)0-2 gf.^^ \V1V2V3V4, 
By performing the derivative with respect to A;2,;ti we can rewrite the integrand as follows 
d f ku \ dku 1 



dk2„ 



V^V2V-iV4, 

^ '-(^]- 

VlVi dk2.t, \V2V3J 
2k ■k2 2k- ki 



1 



2k ■ k2 



ViVlVsVi ViV2V?^Vi ViV2VjV4 ■ 

After expressing the scalar product k ■ ki in terms of the propagators Di 



' and V4 = (ki - k)"^ 



2k ■ ki 

and substituting Eq. II20II in Eq. I|19|l . we obtain, 



■-Vi -V4 + k^, 



I 



d^ki d°k2 



1 



(27r)0-2 (27r)0-2 DiDaCf 



-[- 

J (27V 



d^ki d^k2 



(27r)0-2 (27r)0-2 Vi^lV^ 



, r d"k 

" J (27r)D 
'J (2n 



d°ki d°k2 



(27r)0-2 (27r)0-2 ViV2V^V4 
d°ki d^k2 k ■ k2 



- 2 



J (27T 



d^ki d°fc2 



k ■ k2 



(27r)-D-2 (27r)0-2 ViV2V?J)4. 



Eq. l|21|l can be pictorially represented by 



(18) 



(19) 



(20) 



(21) 



(22) 
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where a dot on a propagator line means that the propagator is squared; irreducible scalar products in the numerator 
are explicitly indicated. 



3.4.2. Example 2: IBP-id's 

Let us consider the case of the two-loop three-leg topology of Fig|2] 




P2 - k: 



Fig. 4. Example of a two-loop three-leg four-denominators topology: p| = p| = —m , and (pi + P2) = Q = ~^ 



The integrals associated to this topology can have three irreducible scalar products, arbitrary chosen to be (pi -fci), 
(P2 ■ f^i) and (fci ■ k2). In doing so, Eq. (|8]l, for generic values of the indices rui and n^, reads 



r d"k 

J ~{2n)0 



ki d°k2 d 



(pi-fcl)"l(p2-fcl)"^(fcl-fc2)"3 



(23) 



where Di = + m^, D2 = k'2 , ©a = (p2 — ^2)^ + m^, 'D4 = (pi + P2 — k-i — k2)'^ + and « = 1, 2. By setting, for 
simplicity, mi = • ■ ■ = m4 = 1, ni = ■ ■ • = n3 = 0, t)p = pi and i = 1, Eq. I I23I I becomes 



r d^ki dOk2 d J- pi,^ \ ^ 



(24) 



After taking the derivative with respect to ki and simplifying the reducible scalar products against the corresponding 
propagators, one can write Eq. ll24ll as follows 
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3.4.3. Example: Ll-id's 

The Lorentz Invariance identity II12I I for the integral of the topology in Fig. |4] with ni = n2 = na = and mi = m2 
m-j = m4 = 1, reads, 



= 





(26) 



3.4.4. Example: Symmetry relations 

One can exploit the invariance of the integrals belonging to the topology of Fig. U 



/ 



d^fcl dgfc2 (pi ■ A;i)"i (p2 ■ fci)"2 (fci . fc2)"3 

T-lTll 7-)™2 7^»"3 7^'"4 
^1 ^^2 ^3 



(27r)0-2 (27r)D-2 15™! X'™^!)™^©™^ ' ^^^^ 

under the redefinition of the momenta running in the nested electron-loop. In fact the two denominators with mo- 
mentum k\ and (pi + p2 — fci — k2), i.e. ©i = fcj + m? and X'4 = (pi + p2 ^ fci — ^2) + "i^ have the same mass 
Therefore the following redefinition of the integration momentum 



A:i = Pi +P2 - fc'i - k2, 



(28) 



that consists in the interchange of the two denominators in the closed electron loop, does not affect the topology of 
the integral; nevertheless the explicit form of the integrand can change generating non trivial identities. Taking for 
instance ni = n2 = = and m\ = m2 = m4 = 1 and ms = 2 in Eq. I|27|l . the substitution ||28|I gives the following 
very simple relation 





(29) 



By choosing ni = n2 = 0, 713 = 1 and mi = m2 = m4 = 1 and = 2, we get a less trivial identity 





(fcl ■ fc2) 



(P2 ■ '^l) 




(30) 
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3.4.5. Master Integrals 

For each of the N{It^r,s) integrals of the class It.,r,s (see Eq.(j4])) we dispose of IBP-id's, Ll-id's and 
symmetry relations, involving integrals of the families up to /j For r — s ~ the nmnber 

of all the integrals involved in the identities exceeds the number of equations; but, in writing down 
systematically all the equations for increasing values of r and s the number of equations grows faster 
than the number of the integrals, g enera ting an apparently overconstrained system of equations — 
as realized first by Stefano Laportall^. After its generation one is left with the problem of solving 
such a linear system of identities which is trivial in pr incip le, but algebraically very lengthy, so that 
some organization is needed for achieving the solutiorff^. 

To this end, one can order the integrals in an appropriate way. In particular a "weight" is assigned to 
each integral: the weight can be almost any increasing function of the indices rii and rrij , such that 
integrals with higher indices have bigger weights. The system can be solved by Gauss' substitution 
rule, by considering the equations of the system one by one, and using each equation for expressing 
the integral with the highest weight in terms of the integrals of lower weight. Then the result 
is substituted in the leftover equations. The algorithm, by now known as Laporta algorithm, is 
straightforward, but its e xecu tion requires a great amount of algebra, which has been implemented 
in several computer codei^^. 

One finds that several equations are identically satisfied, and the original unknown integrals are 
expressed in terms of very few independent integrals, the master integrals (MPs). In doing so, the 
resulting Mi's correspond to the integrals of lowest weight; but as the choice of the weight is to a 
large extent arbitrary, there is also some freedom in the choice of the integrals to pick up as actual 
Mi's (not in their number of course!). 

There are several cases in which more than one MI is found for a given topology, while sometimes 
only one MI is present. It may also happen that no MI for the considered topology is left over in 
the reduction — i.e. all the integrals corresponding to a given topology, with t propagators, can be 
expressed in terms of the Mi's of its subtopologies with {t — 1), and/or less, propagators. In such a 
case one speaks about reducible topologies. 

Finally, three points are worth mentioning. 
First, any given topology with t propagators has {t — 1) subtopologies with {t — 1) propagators, 
{t — l){t — 2) subtopologies with {t — 2) propagators etc., down to as many propagators as the 
number of loops. It turns out, however, that many subtopologies are in fact equivalent, up to a 
translation of the loop variable, and the subtopologies coming from the contribution of different 
graphs do overlap to a great extent. For these reasons the compilation of the system for the actual 
number of all independent topologies, namely those that cannot be transformed one into the others 
by a redefinition of internal momenta, is relatively small. 

The second remark is about the independence of the identities for a given topology. Under an 
idealistic point of view, it might happen that the infinite set if IBP-id's, obtained by considering the 
infinite choices of indices of denominators and scalar products for the considered topology, plus the 
infinite set of IBP-id's of all the parent topologies containing it as a subtopology, could include Ll-id's 
and symmetry relations as subset. Certainly, at the practical level, since one is working with finite 
sets of indices, Ll-id's and symmetry relations can be considered as additional to the IBP-id's, which 
surely speed up the solution of the system and work as further checks of the reduction procedure. 
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Third, we are not able to prove analytically that the Mi's one finds are really the minimal set of Mi's, 
i.e. that they are independent from each other; in any case the final number of the Mi's is quite small, 
so that reducing the several hundred integrals occurring in a typical calculation to a few of them 
is after all a great progress. Studies on the apriori determinatio n of whether a topology could have 
or not Mi's, have been performed by Baikov and Smirnov [SH. Although a general algorithm 
is still lacking, we think that the analysis of leading and subleading Landau singularities (see 
for an extensive treatment) could be related to this issue: since a reducible topology without Mi's is 
expressible in terms of the Mi's of its subtopologies, it should mean that the leading singularities of 
the reducible topology are not independent of the subleading ones, which are leading singularities 
for the subtopologies. 

As a last remark, let us recall that, at the end of the reduction, there is some freedom for choosing 
the basis of Mi's and usually the choice is in general motivated by convenience. For example the 
behaviour of the functions in the D-to-A expansion might determine to select: i) simpler integrands, 
in view of a successive analytic computation; ii) more complicated integrands, but with a better 
convergence, should their numerical evaluation be of interest. 



4. Differential Equations for Master Integrals 

The outcome of the reduction procedure, previously discussed, is a collection of identities thanks to 
which any expression, demanding originally for the evaluation of a very large number of integrals, is 
simplified and written as linear combination of few Mi's with rational coefficients. The completion 
of the analytic achievement of the result proceeds with the evaluation of the yet unknown Mi's. As 
we will see in a moment, the same collection of identities is as well necessary to write Differential 
Equations satisfied by the Mi's. 



4.1. System of differential equations 

Once all the Mi's of a given topology arc identified, the problem of their calculation arises. Exactly at 
this stage of the computation, differential equations enter the game. The use of differential equations 
in one of the internal masses was first proposed out by KotikoA^, then extended to more general 
differential equations in any of Mandelstam variables by RemiddiH 

Let us point out the basic idea of the method. To begin with, consider any scalar integral defined as 

M{s,,s,r--,su)^ j (31) 

where {si, S2, • ■ ■ , s^} is any set of kinematic invariants of the topology and M is the number of 
such invariants. 

Let us denote the set {si, S2, • • ■ , Sf^} = s and consider the following quantities 

0,fc(s)=p,,^^^ (i,A: = l,2,... ,5-1), (32) 

where g — 1 is the number of independent external momenta. By the chain differentiation rule we 
have 
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^ dsa dM{s) 



^ f d.sa,\ dM(s) 



dSa 



(33) 



According to the available number of the kinematic invariants, the r.h.s. of Eq. (|32p and the r.h.s. 
of Eq. ([33)1 may be equated to form the following system 



AT 



9sa N dM{s) 



dM{s) 

dp. 



(34) 



which can be solved in order to express ^^J'^^ in terms of pj^f^ ^dpk^^ ' ^'^^ corresponding identity, 
can be finally read as a differential equation for M. 
Examples of such equations are the following. 

• 2-point case. 

• Differentiation with respect to a mass 



d 
dm?- 



p — — p 



(35) 



where, for simplicity, we assumed there is just one propagator of mass m. 
• Differentiation with respect to the squared momentum 



2_d_ 

dp^ 



1 d 



(36) 



• 3-point case. 
p. 9 



Qp2 



Pi 



P2 



PS 



A, 9 d 

A pi.p^ +P2,M^ 



B[pi^^^^+P2,^. 



with P — pi + P2 and A, B rational coefficients. 
• 4-point case. 



Pi 



P2 



P3 



(37) 



P 



d 



gp2 



















^ P2 









d 



d 



dpi,f, 
d 



P3 



d 



dp3,i^ dpi.^ 9p2„ 



Pl P3 



Dp: 



dp: 



P2 



P4 



(38) 
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,P3. 






M. 




H 


^ P2 









d d d 



d d \ ^ d 
T. P2,m;^ + <^P2,/iT; 

OPl.p. Op2.j,J Op2,f, 

pi P3, 

(39) 



3,M 



P2 Pi 



with P — pi + p2, Q ^ Pi — Pi and C, Z?, F, G, H rational coefficients. 

Equation ([M)) holds for any function M (s). In particular, let us assume that M(s) is a master integral. 
We can now substitute the expression of M in the r.h.s. of one of the Eqs. p6ll55)l . according to the 
case, and perform the direct differentiation of the integrand. It is clear that we obtain a combination 
of several integrals, all belonging to the same topology as M . Therefore, we can use the solutions of 
the IBP-id's, Ll-id's and other identities for that topology and express everything in terms of the 
Mi's of the considered topology (and its subtopologies). If there is more than one MI, the procedure 
can be repeated for all of them as well. In so doing, one obtains a system of linear differential 
equations in s for M and for the other Mi's (if any), expressing their s-derivatives in terms of the 
Mi's of the considered topology and of its subtopologies. 

The system is formed by a set of first-order differential equations (ODE), one for each MI, say 
Mj , whose general structure reads like the following, 

^M,{D,s) = Y,Ak{D,s) Mk{D,s) + J2BhiD,s) iV,,(i?,s) (40) 

where a = 1, • • • ,J\f, is the label of the invariants, and Nk are Mi's of the subtopologies. Note 
that the above equations are exact in D, and the coefficients Ak,Bk are rational factors whose 
singularities represent the thresholds and the pseudothresholds of the solution. 
The system of equations pO)) for AIj is not homogeneous, as they may involve Mi's of subtopologies. 
It is therefore natural to proceed bottom-up, starting from the equations for the Mi's of the simplest 
topologies (i.e. with less denominators), solving those equations and using the results within the 
equations for the Mi's of the more complicated topologies with additional propagators, whose non- 
homogeneous part can then be considered as known. 



4.2. Boundary conditions 

The coefficients of the differential equations (|40ll are in general singular at some kinematic points 
(thresholds and pseudothresholds), and correspondingly, the solutions of the equations can show 
singular behaviours in those points, while the unknown integral might have not. The boundary 
conditions for the differential equations are found by exploiting the known analytical properties of 
the Mi's under consideration, imposing the regularity or the finiteness of the solution at the pseudo- 
thresholds of the MI. This qualitative information is sufficient for the quantitative determination 
of the otherwise arbitrary integration constants, which naturally arise when solving a system of 
differential equations. 
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4.3. Laurent expansion around D — 4 

The system of differential equations (|40|) is exact in D, and, when possible, its solutions could be 
found for arbitrary value of the dimensional parameter. Very often, the result of the computation 
will have to be anyway expanded around some typical value of the space-time dimension. Indeed, 
in what follows, we will discuss the Laurent expansion of the solutions around _D = 4, though the 
procedure can be applied equivalently for other values of D. In general, we use the ansatz 

n 

Mj{D,s)^ J2 (^-4)'= Mf'(s)+o((i?-4)("+i)) , (41) 

fc— — no 

where uq (positive) corresponds to the highest pole, and n is the required order in {D — 4). 

When expanding systematically in {D — 4) all the Mi's (including those appearing in the non- 
homogeneous part) and all the D-dependent terms of (|4T|) . one obtains a system of chained equations 
for the Laurent coefficients ikfj'^'' of ((4T|) . The first equation, corresponding to the highest pole involves 
only the coefficient Alj^ ""^ as unknown; the next equation, corresponding to the next pole in (D — 4), 
involves the Mj^ "o+i) unknown, but it may involve in the non-homogeneous term (if it 

appears as multiplied by any power of (D — 4)); but such a term, AI^ has to be considered known 
once the equation for the highest pole has been solved. For the subsequent equations we have the 
same structure: at a given order fc in (£) — 4), the equations involve Mj'^' as unknown, and previous 

coefficients Mj^-* (—no < i < k) as known non-homogeneous terms. 

Let us note that the homogeneous part of all the equations arising from the D — > 4 expansion of 
(|1T|) is always the same and obviously identical to the homogeneous part of Eq. fiO]) read at D = 4. It 
is, therefore, natural to look for the solutions of the chained non-homogeneous equations by means of 
Euler's method of the variation of the constants, using repeatedly the solutions of the homogeneous 
equation as integration kernel, as we will see in the following chapters. 

General algorithms for the solution of the homogeneous equations are not available; it turns 
out however that in all the considered cases the homogeneous equations at D — A have (almost 
trivial) solutions, so that Euler's formula can immediately be written. With suitable changes of 
variable, according to the kinematics of the problem, all integrations can further be carried out in 
closed analytic form, by exploiting the shuffle algebra induced by integration-by-parts on the nested 
integral representation of the solution. 



5. Laplace Transform of Difference Equations 

Difference Equations are functional relations among values of functions shifted by integers, and 
can be considered a "discretization" of differential equations. The set of identities used to derive 
Differential Equations in the external invariants of Feynman integral, can be used as well to derive 
Difference Equation in one of the denominator powers, as discussed by 

LaportjPIMI. 

Alternatively, 

the shifted dimension relations among scalar integrals of TarasoJ^SI naturally yield equations where 
the integer variable is the dimensional parameter D. We briefly describe the former kind of difference 
equation and their link to differential equations. 
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Let us consider the difference equation 

N 

J2pkin)U{n + k)^0 , (42) 

k=Q 

where A'' is the equation order, and Pk{n) are polynomial in n of maximum degree P, whose generic 
structure can be parametrised as follows, 

P i-l 

Pk{n)^AkQ+J2^kt'[[{n + k + j). (43) 

i=l j=0 

The Laplace transform method consists in the substitution 



C/(n) = / dt f"' v{t) , (44) 

where 7 is a suitable integration path, and whose effect can be shown to be the translation of 
the difference equation p2)) . into a differential equation for v{t), 

E'f.w i-ty ^"W = o, (45) 

4=0 

with 

JV 

<I>,(t) = EAfe, . (46) 

k=0 

The solution v{t) of Eq. ([l5l) . still carry unknown integration constants which will be present as well 
in U(ri), reconstructed by integrating Eq. (|44|) . Finally, the value of the yet undetermined integration 
constants can be fixed by imposing boundary conditions in the large-n regime of U(n), easily found 
by direct inspection of its original representation as loop integral^. 

71 

5.0.1. Example: one-loop Tadpole ..^""^ 

t/(„) = .-f / f\ = LJ . (47) 

By means IBP Id's, one can write the following relationship between the integrals of family II47I I 

71—1 n 



- (^n - 1 - j +(«-!) = 0, (48) 

which, after renaming n — > n + 1, can be put in the form I I42I I. and read as a first-order difference equation for the 
function U{n), 

- (n - — ]U{n) +n U{n + 1) = 0. (49) 



In the rest of the paper, for computational convenience, we adopt two different definitions of the loop measure: 
(i) k / ~^ , for integrals treated with the Differential Equations method; (ii) d^ k/ {tt)^ /'^ , in case of use of 
Difference Equations method. 
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The boundary condition is determined by the asymptotic behaviour (n — » oo) of the integral I I47I I. which reads 

n-- / — —X7T-- / d^ke-"'' =n--. 50 

J (fc2 + l)" J ^ ' 

One way to solve equation I l49t is to seek for the solution as a, factorial series (see | 65 | 66 | Jqj. details) . Equivalently, 
one can convert the difference equation into a differential equation, by means of the Laplace transform of U{n). 
Accordingly, with the ansatz 

U(n) = f dt t"-^ v{t) , (51) 
Jo 

and 

A D 

Am = — , Aoi = -1 , 

Am = -1 , All = 1 , (52) 
equation I I49II becomes a differential equation in the form of (145 l l for the function v{t), 

-t{t-l)v'{t)+(^^^-ty{t)=0. (53) 

The solution is, 

v(t)=vot-^/^ il-t)^/^-'- , (54) 
where vq is the yet unknown integration constant. With the above solution, Ea JSll l can be integrated and U{n) reads, 

U{n)=vo ^ \\ ^. (55) 

r(n) 

By comparing its large-n expansion, 

U(n) "=°°„n n--D/2fr/^ ^ 
with the asymptotic limit explicitly computed in Eq. llSOI I. one finally determines the value of the constant. 



Uin)''=°°von-^/^(r(^)+o(-]] (56) 



6. One-Loop Vacuum Polarization in QED 




p — k 

Fig. 5. 1-loop vacuum polarization diagram. 

As an application, we calculate the one-loop correction to the photon propagator: the so called 
Vacuum Polarization tensor. The only contributing diagram is shown in Fig. [5l and corresponds to 
the expression 

i(-)u (k)-( zc)^( 1) f ^""P Trf ^^^''^+"'^^'^[~'^^~^^"^"'n (581 
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n^i^ is a tensor which depends only on the external momentum and the metric tensor S^i, and 
can be decomposed as a sum of two contributions 

U^,{k) = Il{D, k^)k^k, + AiD, k^)S^,, (59) 

where n(D, fc^) and A(_D, k'^) are scalar functions. The Ward identity (current conservation) tells us 
that kf^Ilf^i,{k) = 0, hence 

k^,Uf,^{k) = =^ k^k'^niD, k^) + k^A{D, fc2) ^ =^ A{D, k^) = -k^Il{D, k^), (60) 
and we can rewrite Eq. (I59p as 



n^,(fc) = n(i?, k^)ik^k, - k's^,) . (61) 

Let us extract the scalar function n(D, fc^) taking the trace of the (pT|) 



n^^i?, fc) = u{D, k'){k' - k^D) =^ n{D, fc2) = __L_n^^(fc). (62) 



. . 2 

From (|58|) . setting a = |-, we have 



(27r)-0-2 (p2 _|_ „j2-)[(-p _ _|_ ^: 
after some operations within the Dirac algebra, our expression becomes 



= -Paqi3"^T:{lt,lalvli}) + "i^Tr(7^7y) = 

= -Tr(I)D Paqf3iSfj.aSvl3 - 5^i.i,5al3 + '5p/3'5ya) + Tr(I)r, 5 
= -Tr(I)_D {PiiQiy - P ■ q5fj.,y + Pi,qi_i - TO^5^i/) = 

= -Tr(I)_D [2p^p^ - Pf^k^ - p^kf^ - {p^ - p ■ k + m^)8^^], (64) 

where Tr(I)D is the trace of the unit matrix in D dimensions (this is, in general, different from D). 
However, lim£)^4 Tr(I)£) — 4. Substituting (I64|) in (I63p and putting fi = v we obtain 



d'^p 2p2 - 2p • fc - D(p2 -p-k + rn^ 



= Tr(I)i3 ^ 

- 1T[1)D I I ^o_\_D-2 (-^2 _i_ ^2Arr„ _ )L\2 _i_ „-,21 ' ^^^^ 



(27r)^-2 (p2 +to2)[(p- fc)2 +m2] 

rf^p (i:>-2)(p-fc)-(i:»- 2)^2-1)1 

(27r)^-2 (p2 _^ j^2)[(p _ /,)2 _^ „^2 



We can define, 



2?1 p2 _|_ j^-, 2 

I?2 = (j3 - fc)^ + 
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SO that after the trivial tensor reduction, according to 

2 7-i2 2 

p = Di — m 



Eq. (|65|) becomes 



(66) 
(67) 



Tr(I)c z 



d^p {D - 2){p ■ k) - {D - 2)p2 _ D 1 



Tt(1)d i 



Tr(I)z5 i 



d^p [D - 2)[(X>i - X>2 + fc^)/2] - (£> - 2)(Pi - m2) - Dm^ 

(2-1?) /- d^p f I J_\ , [fc2(L»-2) 
(2^)^-2 Iv^ + P^j ' 
P(i:)-2) 



2 

(2-1) 



2to" 



d^p 1 
(2^)^-2 p^p^ 



■(Jio + Joi) + 



2777, 



where all the integrals can be seen as belonging to the class, 

d^p 1 



(27r)^-2 P^ip™^ ■ 

Finally, the scalar vacuum polarization function can be written as 

n{D,k'^) = 



2A;2(1 -D) 
1^(1)^ I 



{(2 - i?)(Jio + Joi) + [k\D - 2) - 4777^] Jn} = 
_ ^) I (2 - D)T{D,m') + ^[k'{D - 2) - 4777^] J(Z?, /c^^ 



One may notice the result is expressed in terms of two integrals: the massive Tadpole 



m 



T{D, m^) = Jio = Joi = 
and the 1-loop 2-point function 

J{D,k^) = Jn 

where 



d^p 1 



C{D)- 



,D-2 



{2tt)D-2Vi ' '{D-2){D-A) ' 
d^p 1 



C{D) = (47r) — r(3-i:»/2) 



(68) 



(69) 



(70) 



(71) 



(72) 



(73) 



is an overall loop factor whose value in 4 dimension is C(4) = 1. In the following pages, when not 
ambiguous, we will give as understood a factor C{D) for every loop. 
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6.1. Reduction to master integrals 

We have seen that the scalar vacuum polarization function Il{D, k"^) is reducible to the sum of a 
Tadpole and an integral belonging to the 1-loop 2-point topology, the most general integral of this 
class being 



d^p 



1 



By means of IBP id's, we can reduce every integral of the type ([7i)l to 

d^p 1 



(74) 



(75) 



which is, actually, the only MI of this topology, and to the tadpole Jio, which is the MI of the only 
possible subtopology. 

Let us write explicitly some useful IBP id's for the integral (|75p. we have 



r ° P 



d"p d 







(76) 



where I>i = + and I?2 = {p — k)^ 
identities 




with 



(D-3) 
'(4m^+/c^) 



(4m^+fc^) 



(77) 





1 


(27r)^- 






1 


(27r)^- 






1 


(27r)^- 




d^p 


1 


(27r)^- 





(78) 



(79) 



(80) 



(81) 
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where the last equation has been obtained as IBP identity for the tadpole 



m 



6.2. Differential equation for J{D,k^) 

The master integral J{D, k^) is an analytic function of the argument fc^ and it can be viewed as the 
solution of a suitable differential equation. Let us see how to build and solve such an equation. For 
J{D,k^) the following trivial identity holds, 



dJ 

dku. 



dJ dP , dJ 



(82) 



By contracting 



with the vector fc^ we have 



On the other hand 



, dJ__ „,2 9J 



(83) 



dJ 
dku 



d^p d 



{2n 



i_D-2 



dkf, \V1V2 



{2n 



(84) 




By substituting Eq. 



in Eq. 



we have 



(85) 



(86) 



which is rewritten, thanks to the second identity of the ([77|l and to ([8T|) . as a non-homogeneous 
first-order differential equation for J{D, k"^) 



1 (i? - 3) 



(£-2) 
Am? 



k"^ {kP-+Am?) 
1 



fc2 (P + 4to2) 



(87) 



Eq. (j87p contains the boundary condition for the solution. In fact, thanks to the analytic properties 
of Feynman integrals, we know that J(-D, k"^) must be a regular function in fc^ ~ 0, that is 
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hm fc2 — = 0. 

By multiplying Eq. ([87|l by fc^ and taking the limit fc^ — > 0, we have 



lim k'^-rriT 



lim — 

fc2^o 2 



1 - 



(L>-3)fc 



2 1 



+ 4to2) 



lim 



(fc2 + 4m2 



out of which one has, 



\im^JiD,k^) = J{D,0) 



(^-2) x 



(89) 



(90) 



Eq. ((87)) with the condition (|90p constitutes the Cauchy problem we have to solve. 



6.3. Exact solution in D dimensions 



Equation (|87p can be solved exactly in D dimensions. Let us introduce the dimensionless variable 
X = obviously — -^-j^ = 4^^i ^^^d the equation ([57)1 can be rewritten as follows 



1 _ (£-3) 
X {l + x) 



{D-2) 
Am? 



X {l + x) 



(91) 



The general solution of Eq. (|9ip is the sum of solution of the homogeneous equation, say Jq, and 
particular solution of the complete equation, say J*. The homogeneous equation is 



dx 



1 (D-?,) 



X (l + x) 



Jo 



with solution, 



Jn{D,x) ^ Ax 2 (1 + x)' 



(92) 



(93) 



To find a particular solution, J*{D,x) of the complete equation, we use Euler's variation of 
constants method, and write 



J*{D,x) =a;-^(l + x)^^0(x) 



(94) 



where (l){x) is an unknown function to be found by imposing J* fulfills Ea. (j9ip . In so doing, we 
obtain the following equation for (^{x), 



dx 



(D_2) 
Am? 



i(l+x) 2 _ a;2 (1 -(- a;)- 



(95) 



hence 
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Am? 



X 2(1 + 2 dx — x-^{l + x) 2 dx 



The integrals in Eq. (|96p are a representation of hypergeometric functions 



hence 



/ x-Hl + x)-^dx^2xi 

'3 D-1 . 5. 

2 ' -^Z ' 



(96) 



(97) 



7(1), x) = Jo + J* = 
= Ax 5(1 + 2;)^— ^ ' 



Am? 



(1 + x)~ X 



„ ^ , 1 D~i 3 



2 ^ /3 5 
—2:2^1 — , \—\~x 

3 V2' 2 '2' 



(98) 



By imposing the boundary condition, we see that the constant ^4 = 0, because the term x~2 (1 - 
a;) * 2 ' is singular for x — *■ 0, while the MI is regular by inspection. 
Finally, the full solution of our Cauchy problem is 



J{D,x) 



(D-2) 



2m^ 
, 1 3 



T 2 '2 
P-2) 



2to^ 
(£-2) 
2m^ 



(1 + x) 



(1 + x)~-x 

X /3 D-l 5 

1^-3 /£>- 1 1 3 
^ ^^U^'2'2'-" 



4-D 



2 '^'2' 



(99) 



6.4. Renormalized vacuum polarization function 

Now we can write the exact i'-dimensional expression for the renormalized vacuum polarization 
function. This is 



n{D,e)R^ii{D,e)-niD,o), 

where Tl{D, 0) is the function evaluated at zero momentum transfer 



(100) 



U{D,0) 
We see that 



Tr(I)D t{2-D) 
(1- D) Am^ 



liml + ^J(AO)-ilim:^. (101) 

1:^0 X 2 ^ ' ^ 2 x^o X ^ ' 
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lim ii£l^ ^ 0) Km - + J'{D, 0) = O 1™ - + 0) , (102) 



hence 



Finally, one has the renormalization counterterm, 



11(0,0) =Tr{l)n^^-^J^ Q (104) 



yielding the renormalised form factor, 



IlRiD,k^)^U{D,k^)-n{D,0) 



X < - 

X 



{D-2)-- 

X 



Tr(I)g i(2-D) 
(1 - D) Am? 

2 ' '2' 



,'4-1? 3 \ 2, 
2F1 ( ^^,l;-;-x +-(l-i?) 



(105) 



which can be ultimately expanded around Z? = 4 up to the finite order in terms of HPL's 



9 3(l-z)' 3(l-z) 
(1 + z) (l-4z + z2) H{0,zy 



3{l-zr 



with the spacelike variable z defined as 



\/—x + 1 — y^—x V— fc^ + 4m^ — V— fc^ 
V-a: + 1 + V-fc^ + 4m2 + V-fc^ 



(106) 



(107) 



7. Two-Loop Vacuum Polarization in QED 

The two-l oop c orrections to the vacuum polarization function was first calculated in 1955 by Kallen 
and Sabr y i At that time, it was one of the first applications of perturbative QED, soon after 
the theory was definitely laid out. The exact e xpre ssion in dimensional regularization was calculated 
in 1993 by Broadhurst, Fleischer and Taraso v I The result was expressed in terms of generalized 
Hypergeometric functions 2F1 (4^) and 3F2 

( 4m^ ) ' ^^*3re k is the momentum of the external 

photon and m the mass of the internal fermion. 

In the following sections we will derive the same result by means of the method of differential 
equations. First of all, we will show that the use of IBP allows one to reduce the whole problem to 
five Mi's, out of which three are product of one-loop integrals (known from the previous section), and 
two are actual two-loop integrals, as yet unknown. Then, we will write for them a system of two first 
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order differential equations in the square of the external momentum. We will see that the system 
is equivalent to a pair of second order differential equations, one for each MI. These equations 
can be solved exactly in D dimensions and the solution can be written in terms of generalized 
Hypergeometric functions. 



7.1. Diagrams 

The one particle irreducible contributions to the two-loop vacuum polarization are showed in Eq. 
pOSp . They consist of three bare diagrams, plus the corresponding counterterms, namely two fermion 
self-energy counterterms and two vertex corrections counterterms. Due to Ward identities of QED, 
some of these counterterms cancel each other: the vertex subtractions cancel exactly the wavefunc- 
tion renormalization counterterms; so the only effective contribution comes from the fermion mass 
counterterm, as shown in Eq. (|108p . The amplitude reads 



i[-^ n^,(fc) 



(108) 



iSm 



7.2. Topologies and Master Integrals 

From the above Feynman diagrams, one can identify seven independent topologies, which cannot be 
related to each other by a transformation of the internal momenta, as shown in the second column 
of Tab. [11 



By the systematic application of the reduction algorithm, one can express all the needed integrals 
as combination of just five Mi's, depicted in the last column of Tab.[TJ one MI with four denominators; 
three with three; and one with two. 



denominators 


Independent topologies 


MI 


5 




completely reducible 


4 






3 






2 


m 


CX) 



By giving a close look at them, one soon realizes that three of them are indeed product of the 
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one-loop Mi's we met in the previous section, namely, 



QQ=T2(Am^), (109) 
--{^y^ = J{D,k^)xT{D,m^) , (110) 
™(^^^^^^(^^^ =J^{D,k^), (111) 



where T{D,m?) and J{D,k'^) were given respectively in Eqs. (TfTl [99l) . 
The two yet unknown Mi's, both belonging to the same topology, are, 



M,') = = j (3^^ (112) 



where Di = p\+ , 1^2 = P^, = (k — pi ~ P2)^ + tti^ , Pi and p2 being the two loop momenta and 
k being the external momentum. Let's discuss their evaluation. 



d^pi f d'^p2 {pi-P2) 



7.3. Differential equations for Ji{k^) and J-2{k^). 

The functions Ji(fc^) and J2{k^) are found to fulfill a system of first-order differential equations in 
the external momentum squared that reads as 



( Ml. — 



iD-3) , (31^-7) 
2fc2 "1" 2(fe2+4m2) 



Jl 



3(D-2) 



j,2 k'^+4n 



■h- 



(g-2) 



1 



^,2 k'^+im^ 



dJ2 _ (g-2) T {D-2) T 

4 'Jl oi,2 'J', 



^ , _ {D-2) rp2 

2fe2 ■J2 4fe2 -l- 



dk^ 

where T{D,m?) is the massive tadpole. 

Such a system is equivalent to a second-order differential for Ji 



77^ (?7 -I- 4771-^) 



2^ ^ZX.- 



(D-3) 



(^-4) 



(Z) - 4) - 6m^ 



drj 



V- 



- {D- 2)m" 



y2 



0, 



(114) 



(115) 



where rj = k^. 

Differentiating once more with respect to 77, we obtain a third order differential equation for Ji 



r/2(ry + 4m2)^ 
(i?-4)(i?-9) 



(I? - 6) 77 - 14m^ 



77 - D{D - 5)™^ 



C?77 



{i-D){A-D) 



Jl =0, 



(116) 



which, introducing the dimensionless variable x 



becomes 
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+ 



(PJi r3(i:»-6) 7 



(L>-4)(7:>-9) D(D-5) 
X ■ 



dJi (3-D)(4-D) 
dec 2 



Ji = 0. 



7.4. Exact solution in D dimensions. 

Quite in general, a differential equation like 



x\l + x)-j^ + x[x{ai + aa + as + 3) + (5i + 62 + l)]^:^ 



dx'^ 



.dip 



+ [x{aia2 + aifla + aafla + ai + 02 + 03 + 1) + ^1^2]-^ + 
+ {aia2a3)ip = 0, 



(117) 



(118) 



where 01,02,03,61,62 are parameters, is classified as hypergeometric and its solution space is 
spanned by the functions 



3^2 



oi, 02, 03 ^ 
61 , 62 



ijji{x) 



1 + Ol — 61, 1 + 02 — 61, 1 + 03 — 61 

2-61,1 + 62-61 



(119) 



1 + Oi - 62, 1 + 02 - 62, 1 + 03 - 62 , _ 
2-62,1 + 61-62 

We can easily see that Eq. pi7p is hypergeometric if we identify the parameters as follows: 



f Oi = 1 

02 = 3 - L» 
03-(4-i?)/2 
61 = D/2 
I 62 = {5-D)/2. 

With the above choice, the three independent solutions of Eq. (|117p become 



D/2,{5~D)/2 ' 

(-P-2) 



ijji{x) = 


3-P2 




i~x) 


^^3(2;) = 


i-x) 



3-F2 



3^2 



B-3I?)/2,(3-i^),(4-i^)/2 _ 
(4-I?)/2,(7-2I?)/2 ' ^ 
(3-Z?),l/2,(i?-l)/2 _ 
(2i:>-3)/2, (i:>-l)/2 ' 



(120) 



(121) 
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Finally, the general solution of Eq. (I117P can be written as a linear combination of Tpi (x) , 7/12 (x) 
e tlj3{x) 



Ji{x) = Aiji{x) + Bi}j2{x) + Ci^zix) 



(122) 



Now we can choose the value of the three integration constants A, B and C by means of suitable 
boundary conditions. By inspection we see that for a; — > the MI Ji becomes 



Ji(0) 



(27r)^-2 (pI + m'^)pl[{pi - P2Y + m'^) 



(27r)^-2 

CO' 

where the last line has been obtained by means of IBP-id's. If we study the behaviour of (I123|) 
in the UV and IR limits, we see that it is regular for 2 < D < 4. This means that B = C = 0. 
If this was not the case, in fact, terms such as {—x) 2 and (— x) 2 would give rise, in the 
abovementioned range of D, to a divergent behaviour, and this would not be compatible with the 
finite resuh of 

Furthermore, (|123p allows one to choose once and for all the constant A 

Ji{Q) = AiPi{0) = A, (124) 
because, like all Hypergeometric functions, ipi{0) = 1. Comparing (|124p with (|123p . we have 



A = 



hence 



Jiix) 



(D-2) 



3^2 



(3-D),(4-D)/2,l 
D/2,{5-D)/2 ' 



(125) 



2to2(d _ 3) 

Now we can obtain the expression for J2{x) by simply substituting Ji{x) and in the first 
equation of (|114p . To do so, we have to know the expression for the first order derivative of the 
Hypergeometric function 3-F2(a^)- 
Quite in general, a function like 



3^2 



has the following series representation 



61,^2 



3-P2 



01,02,03 
61, 62 



(ai)„(o2)„(a3)„ 



(126) 
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where {Cjn is the Pochhammer symbol, defined as 

r(e + n) 



(Or. 



m 



By differentiating Eq. (I126P with respect to x, we obtain 



dx 



3F2 



E 



61,62 

(ai)n(a2)n(a3)n (-2^)'' 



'^^ ^0 (^l)n(^2)n 



E"*""" (ai)K(a2)«(a3)K ^ 
(6i)„(62)„ " nl 



_ _ (Qi)K(a2)n(a3)n (-a:)" ^ _ 
t{ ibiUb2)n in-iy. ~ 

_ (ai)fc+i(a2)fe+i(a3)fe+i {-x)'^ 



fc=0 

Now, we see that 



ibi)k+i{b2)k+i 



k\ 



(127) 



hence 



r(g + fc + i) _ r(e + i + fc) _ 
(?)fc+i = — T^TT^ — ~ ^~FrrT~r\ — ~ + ^i^^ 



r(e + i) 



3-P2 



ai, 02, 03 . 

61,62 ' 



(aiKoaKos) (ai + l)fc(a2 + l)fe(a3 + l)k {-xf 
(6i)(62) 



■E 

fe=0 



(aiKa^Kos) 
(6i)(62) 



■ 3^2 



(6i + l)fe(62 + l)fc 
ai + 1, 02 + 1, as + 1 



k\ 



61 + 1,62 + 1 



Solving pi4p with respect to J2 we have 



(128) 



J2{x) 



i{D - 2) 
2m2 



[2{2D - b)x + [D - ?,)]Ji{x) 



3(D-2) 

Substituting (|125p in (|129p and because of 



a;(x + 1) 



dJi(x) 1 
da; 3 



T'^{D,m'^). 



(129) 



3— 3^2 

ax 



(3-i?),(4-i?)/2,l 
I?/2, (5 - i^)/2 ■ 



2{D - 3){D ^ A) 
D{D-b) 



■ 3F2 



{A-D)A6-D)/2,2^ 
{D + 2)/2,{7-D)/2' 



, (130) 



we can write down finally the exact expression in D for J2{x) 
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CO 



3F2 



(3-2?),(4-2?)/2,l, 



D/2,i5-D)/2 

(4-Z?),(6-Z?)/2,2_^ 
{D + 2)/2,{7~D)/2' ' 



(131) 



7.5. Renormalized vacuum polarization function 

By means of the integrals' reduction, the two-loop Tl2Fi., with the one-loop subdivergences already 
subtracted off according to Ea. (|108p . admits the following decomposition in terms of MPs, 



with 



n2L,lfl(i?, fc ) = Cl -\- C2 -|- C3 

+C4.^/^Y^ -hC5 00 , (132) 



—2 (—2 + L*) 

" (-12+19Z?-8i?2 + ^3)p^2 (A:2+4m2) ( ^"^^ + " + 



-t-4 (-14 + 22D - 9D^ + D^) k^m-" + 16 (-3 + Df j ; (133) 

12{-2 + Df {{8- 5D + D^) P + 4: {10-7 D + D^) m^) 
" (-12 + 19D-8i:'2 + i:»3) fc2m2 (fc2 + 4^2) ; (134) 

+A (-4 + £))^ (-1 + D) - 32 (-3 + D) m"^ ; (135) 

C4 = 71 TT^ ( -2 (-32 + 30D-9DOD^) 

{A- 5 D + D^) k^ (k^ + Am^)\ ^ ' 

-8 (-40 + 38i:i-llL>2 + i:>''^) /c2^2_^64to'') ; (136) 



_ 6(-2 + L>)^ ((8-5L' + L'2) fc2 + 4 (10-7L» + L'2) 

''^ " (-12 +19 £'-8 1)2 +1)3) fc2^2 (fc2+4^2) (137) 

and where we used the mass countcrterm defined, in terms of the one-loop tadpole, as, 

(138) 

2 m (Z? - 3) ' 
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The completion of the two-loop renormalization procedure requires the subtraction of the value 
of n(2L,i_R) at zero momentum, which can be obtained from the above expression, 

n rnn^^ M- d + {-12 + d) d) 

The two-loop renormalized expression 

T12r{D, e) - Hai.iflp, fc2) _ II^l,ir{D, 0) (140) 
agrees with the result in literature I -^^^^ i 

8. System of Three Differential Equations 

After having shown in its completeness the calculation of the vacuum polarization at 1- and 2-loop 
in QED, we proceed with the discussion of a less trivial case of differential equations whose solutions 
ISSI are required to conipute a set o f integrals which enter the 2-loop QCD corrections to the heavy- 
quark form factors I 84 | 85 | 86 | 87 | 106 | g^j.^ going to describe directly the evaluation of three master 
integrals belonging to the same topology. In this case we give as understood i) the reduction to 
the master integrals, namely assuming that all the Mi's have been already identified; and ii) the 
knowledge of the MI belonging to the subtopologies which enters the non-homogeneous term of our 
differential equations. 

We will see, in this example, the technical details related to the Laurent expansion of the system 
of equations around specific values of the dimensional parameter (Z? — > 4), and to the choice of the 
boundary conditions. 

According to our assumptions, we have got the table of identities for reducing all the integrals 
belonging to the topology in FiglHland its subtopologies. 




Fig. 6. Two-loop vertex diagram: = p| = —m?; k = pi + P2', ^1,2 are loop variables; a curly line stands for 
massless propagator; a solid line, for propagator of mass m. 

The scalar integrals represented by the topology in FiglS] are 

J^(ni, 722,723; mi, ma, mg, m4) = y (27r)g-2 V^^V^^V'^'V^^ ' ^ ' 

where Vi = kl+w?,V2 = k^Va = (pi - ki - fc2)^2?4 = {p2 + ki -t- fe)^ 

At the end of the reduction, one is left over with three master integrals, $i(i = 1,2,3), which we 
choose to be 

<^i{D,k^)= (T)-; ^2{D,k^)= (J)^^-' '^3{D,e)= (Pi- . (142) 
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These MI are found to fulfills the following system of first-order ODE, in the variable k^, corre- 
sponding to the momentum transfer, 



d 



(22- 131? + 2 £12) 

{-2 + D) {-7 + 2D) {~8 + 3D) 

i-'S + D) fc2 (fc2 +4^2) 
i-i + D) {-5 + 2D) m2 
2 i-3 + D) (P + 4m2) 
{-4 + D) {-2 + D) (-5-h2D) 

2 (-3 + L») k? {k^+Am?) 
(-4 + 13) {-b + 2D) (-10-f SU) (-8-H3J)) 
4 (-a + Z?) {-1+2D) fc2 (fc2-K4TO2) 



(fci ■ pi) 



-ai) 



(143) 



dfc2 



(fci ■ Pi) (-3 + D) ((-4 + £>) fc2 + 4(-3 + D) m^) 



8fc2 



((20 - 16 £> + 3 D2) p _^ (-5g _ 56 D-t- 12 , 



4fc2 (fc2 -I- 4m2) 



(-4 + 1?) m2 



(fci ■ pi) 



(-236 + 244 D - 82 £12 9 Ll^) 
16 (-7 + 21?) fc2 



(144) 



&2 



(-4 + D) (-3 + D) ((-2 + D) fc2 + 4 (-3 + D) 



4 fc'* m2 

-4 + 1)) {-2 + D) {-8 + 3D) 



(ki -pi) 



2fc4 to2 

((4-6D-h-D2) fc2-|-4 (-2-4i:>-h£)2) ^2 

4fc2 (/c2 + 4 77x2) 

(-2 + 1?) ((14 - 8 D + 1)2) fc2 -I- 4 (8 - 6 L» + D^ 



Ak^m'^ {k^ + Am?) 



3D) 



8 {-7 + 2D) k^m? [k"^ +Am? 
-4 (-104 + 98 £1 - 30 £12 + 3 D^) n/^ 



-164 + 136 £> - 36 132 + 3 D^) P 



(145) 



The non-homogeneous terms of the above equations contain the integrals. 
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J{D,k'^) xT{D,m^) , (146) 
(m2)^-3 T{-D - 5)r2((D - 2)/2)T{2D - 5) 



{D - 4) r(-(D - 6)/2)r{D - 2)r((3L» - 9)/2) ' 



(147) 



which are Mi's of the non- vanishing subtopologies - any other subdiagrams would contain a massless 
tadpole, vanishing in dimensional regularization. The former is the product of the massive tadpole 
T(D,TO^), given in Eq. (ffTj) . and the 1-loop 2-point function J{D,k^), in Eq. ((99|) : while the latter, 
known as 2-loop sunrise, could be easily evaluated by direct parametric integration. 

In principle, solving a system of three first-order ODE's is equivalent to solve a single third-order 
ODE for one of the three MI, say <I>i. But instead of writing directly the third-order ODE, one can 
observe that: in Eas. (|143ll44|) . $3 appears to be multiphed by (D - 4); and in Eq. (|145|) . $1 and $2 
are multiplied by (I? — 4). This features means that, after expanding around D = 4, the original 
system of three coupled equations, get simplified: it decouples in a system of two coupled equations 
for $1 and $2, plus a single equation for $3. 

It is important to remark that the shape of the differential equations depends strongly on the 
choice of the MI: a different choice for $i(i — 1,2,3) could lead to a system which does not get 
simplified after the Laurent expansion either. In all the problems we studied, given the arbitrariness 
of the choice, the practical criterion of having to deal with a simpler system of differential equations 
has determined which integrals had to be picked up as master ones - though there is no apriori 
assurance to find any simplification at all. 

8.1. Laurent expansion in {D — 4) 

Indeed, one rearranges the system of equations (|143ll44ll45p as follows. 

• Writes a second-order ODE for <i>i, from Eqs. (|143ll44p . 

• Take the first-order ODE Eq.(Il45l) for $3. 

• Perform the change of variable 



k^ ^-== --^=^= . (148) 

\/—k'^ + y —k? + rn?- 

Finally, expands both equations around Z? = 4, with a Laurent ansatz for the solutions, 

00 

$1 (7?, x) = ^{D~\f ^^^^ {x) (149) 

fe=-2 

00 

$2(Aa;)= 5](I?-4)'=<i>f (x) (150) 

00 
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After the Laurents expansion, the first-order ODE for <i>3(-D, x) will induce, order-by-order in (Z? — 4), 
a system of chained first-order ODE for the Laurent coeffiecients, <^^j^\x)^ which at the first two 



orders, k - 




'1, reads, 



1 



d 

dx 



= 



(-1 + x {1 + x) 
1 + X iJ(0, x) 



(-2) 



(152) 



H{l.x) 



8 (-1 + x)' 8 (-l + x) (1 + x) 4 (-1 + x) (1 + x) 
(1 + x) (l + x2) $(~2)(x) 4 (1 + x) $^"^^(x) (l + 6x + x2) $^-^)(x) 



2(-l + x) X 
(l + x^) 
(-1 + x) X (1 + x) 



(-1 + x)^ 



2 (-1 + x) X (1 + x) 



d_ 

dx 



(153) 



We remark that the non- homogeneous term of (|153p . the equation for $3 ^''(x), demands for the 
knowledge of the previous order solutions, <i>^ '^\x), <&2 $3^^-'(x). 

Analogously, for the second-order ODE for $i(L),x), the system of chained second-order ODE 
just for the first two coefficients of the Laurent series of <i>i(I?, x), for k = —2, —1 reads, 







1 



= 



4(-l + x) X 
-1 



(-1 + x) X 
H{0,x) 



X dx dx^ 



><S>['^\x) 



(154) 



H{l,x) 



(-2) 



(x) 



8(-l + x)^x 4(-l + x)^x 2(-l + x)^x 



2x2 



(l+x2) (l+6x + x2) d 



(-1 + x) X2 

2 



2 X — 2 x^ dx 



(-2) 



(x) 



(-l + x)^x 1 + x fix dx'^ 



(155) 



In this case, the non-homogeneous term of (|155p . the equation for <i>^ requires the knowledge 



of the previous order solutions, ^[~^\x) and $3~^''(x). 

Therefore, by looking at the Eqs. (ll52|154p for k — —2, and at the structure of the non- 
homogeneous term of Eq. (|153|155|) for fc = — 1, the computational strategy is soon outlined: 

• fc = -2. 

(1) solve Eq.dinil) to find ^i~^\x); 

(2) solve Eq.HMl) to find $["^^(x); 

(3) invert the D ^ A expansion Eq. (|143[) . and substitute the expressions of $3 ^''(x) and 
$i~^^(x) in it, to find ^2~^\x). 



,(-2). 



fc = -1. 



(1) plug the results of the previous order, ^[ $2 '^'i^) ^3 ^''(^)j i'^ the non- 

homogeneous term of Eqs. (|153|155p : 
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(2) solve Eq.HSHl) to find ^\x); 

(3) solve Eq. p55)) to find ^']~^\x); 

(4) invert tlie D ^ A expansion Eq. p43p . and substitute the expressions of $3"^^ (a;) and 



$i \x) in it, to find $^ '{ 
The construction of Laurent coefficients for /c > goes on by repeating the steps (1-4), 
• j, j > 0. 

(1) plug the results of the previous orders, <i>^*''(x), ^2\^) ^-i^d '^'3*''(^) (~2 < i < j), in the 
non-homogeneous term of the equations; 

(2) solve the first order equation for $3"'^(x); 

(3) solve the second order equation for <i>^"''(a;); 

(4) invert Eq.^^ to find $^^^(a;). 

Let us see in detail the case k = —2. 



8.2. Homogeneous differential equations 

The homogeneous differential equation is the same at every order in (_D — 4) , as one can realize by 
looking at the operator in the curly brakets in Eas. (|152ll53p for <i>3*^''(a;) and, correspondingly, at 
Eqs. (|l54ll55|l for $f^(a;). 

The solution 03 (x) of the first-order homogeneous equation is 

X 



(l-x) (l + x) 



The two solutions 0i.i(x), 0i,2(a;) of the 2' order homogeneous equation are 

X 



•^1,2(2;) 



(1 -x) (1 + x) 

X / 1 

(1 - a;) (1 + x) \x 



x-2H{0,x] 



whose Wronskian reads 



W{x) 



(f>l{x) 02 (x) 



-{i + xY 



(156) 

(157) 
(158) 

(159) 



8.3. Solution of the chained equations 

Let us discuss the solutions of the system of equations at order k = —2, according to the above 
strategy. 



(1) One easily solves Eq. (|152p . finding 



l-x 



2^3 



(160) 



(1 —2) 

where $3 is an integration constant to be later fixed by imposing the boundary conditions. 
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(2) Afterwards, one solves Eq. (|154p by Euler variation of constants, 



(2,-2) 



w{x'y 

^ dx' 



Wix') 



(161) 



with <i>^*' = 1, 2) being integration constants, and '^\x) the non-homogeneous term of 

Ea. p54p . thus obtaining. 



^(-2)(^)_ x{l + H{0,x)) 



(1,-2) 



4 (-l + x) (1 + x) (-l + a;)(l + x) 1 
(-l + x^ -2a:g(0,a;)) _^(2,_2) 



{-l + x) [l + x) 



(162) 



(3) Then, the knowledge of $3 '^\x) and (a;), and of theirs derivatives, allows the determination 
of the the coefficient of Laurent expansion of the third MI, by inverting the D ^ A expansion 
of Eq.dHSl), 



^r\x) 



(l + x)(l + g(0,a;)) ^ {l + x) ^(i,_2) 



32 (-l + x) 



(-l + x) 1 



(1+.T) (-l + X^-2.Tg(0,x)) ^(2,^2) 

S{-l + x)x ^ 



(163) 



At the end of this three steps, one knows the expressions of the 1/(_D — 4) ^-coefficient of the three MI, 
'^i{D,x),{i = 1,3) up to the determination of the reaZintegration constants, (f>3 ' , $i , $i ' 



8.4. Boundary conditions 

Usually boundary conditions for Feynman integrals can be read at special values of the kinematic 
variables, and they do correspond to integrals belonging to subdiagrams, therefore to simpler func- 
tions. That is true when the limit to such a specific values is smooth, as it happens around the pseudo- 
thresholds of the corresponding diagrams. In our case the value of x = — 1, meaning — 4m^, is a 
pseudo-threshold. To reach the point x = — 1, being x a space-like variable, as defined in Ea. (|148p . 
one need an analytic continuation to the region, x ^ y — —x + ze, where the Mi's develop an 
imaginary part. 



(-2) 



(-2) 



iv) 



(y) 



y) (i + y) 
-Hio,y)) 



(1,-2) 



(164) 



(1,-2) 



4(-l + y)(l + y) (-l + y)(l + y) 
(-l + 7/ + 2yi/(0,y))_^(2,_2) _y 



(l + 8$^' 



(2,-2) 



i-l + y) (l + y) 



4 (-l + y) (l + y) 



(165) 
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Their expansion around y = 1 reads, 

= 2(l-i,) '^^'''" + 0[(1- vf) (166) 

Finally, the three conditions needed to fix the value of the arbitrary constants, order by order in 
— 4, are: i) the regularity of ^'"^^y) as y ^ 1; ii) the regularity of ^[''\y) as y 1, meaning 
the vanishing of the 1/(1 — ?/)-coe5icients in both cases; Hi) the realness of the constants, meaning 
that real and imaginary part of the 1/(1 — 2/)-coefficient must vanish separately. From the above 
equations, that translates to, 

$^^-^)=0; (168) 

We have all the ingredient to go up in the chain of Laurent coefficients, and repeat the previous 
steps for the case k = —1. The iterative structure of the method yields a bottom- up reconstruction 
of the three master integrals, ^i{D,x), $2(-D,x), and $3(1?, x), around D — A. 



9. System of Four Differential Equations 

This section is devoted to the analytic evaluation of the Mi's associated to the 4-loop sunrise graph 
with two massless lines, two massive lines of equal mass M, another massive line ofmass m, with 
m ^ M, and the external invariant timelike and equal to m^, as depicted in Fig.[7]fl^. 



P 




Fig. 7. Four-loop sunrise diagram: P^ = —rr?\ fci 2,3,4 are loop variables; a wavy line stands for massless propagator; 
a solid line, for propagator of mass M; a dashed line, for propagator of mass m; 



In this case, the heart of the analytic calculation is the study of a homogeneous fourth order 
differential equation, whose solutions turned out to be, in a remarkably simple way, either a rational 
fraction or repeated quadratures of rational fractions. The required four-loop integral could then be 
obtained almost mechanically by repeated quadratures in terms of HPL's. 

Following the reduction algorithm - by now sounding familiar to the reader -, we identify the 
Mi's of the problem; write the system of differential equations in a; = m/M satisfied by the Mi's; 
convert it into a higher order differential equation for a single MI; Laurent-expand it around _D = 4; 
solve the associated homogeneous equation aX D = A and then use recursively Euler's method of the 
variation of the constants for obtaining the coefficient of the [D — 4)-expansion in closed analytic 
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form. The result involves HPL's of argument x and weight increasing with the order in (13 — 4). The 
integration constants are fixed at a; = . After solving the differential equations for arbitrary x, we 
will show how to compute, independ entl y, the numerical value of the solution at a; = 1 by using the 
Finite Difference method of Laporta discussed in SeclHl to show the relation among Differential 
and Difference Equations for Feynman integrals. 



9.1. Master Integrals and Differential Equations 

We find that the problem has 5 Mi's, which we choose to be 

fd°ki...d^ki 

' ~ J (2^)4(C-2) ' ^ ' 

where the numerators Ni are (M'^, ki ■ k^^p ■ fca, fci • k2,p ■ k2) and the denominators I?i = fc^, T>2 — 
I?3 = fc| + A'p, ~ kl + Af ^, P5 — {P — ki — k2 — ks — ki) + .In terms of the dimensionless 
variable x = m/M and putting P ^ Mp one can introduce 5 dimensionless functions ^i{d,x) 
through 

F,{D) = M^^-s C^{D) ^,{D, x) , (170) 

4— D 

where C{D) = (47r)T~r(3 — D/2) is an overall loop normalization factor, with the limiting value 
C(4) = 1 at D = 4. 

The derivation of the system of differential equations is straightforward; the derivatives of the 
Mi's, i.e. of the 5 functions ^i{D,x), with respect to x are easily carried out in their representation 
as loop-integrals Ea. (|169p : when the result is in turn expressed in terms of the same Mi's, one obtains 
the following linear system of first order differential equations in x 

d<^,{D,x) _ j{3D^7) , 3 {D~2) 3 {D ~ 2)\^^^^^^^ _ j 3(0 - 2) , 3 (0-2) 



dx X 

d^3iD,x) _ j3{D-2) 3(0-2) 
dx ^ ^1 2 (1 -x) ^ 2 (l + x) 

j3{D-2) 9(L>-2) 9(L»-2) 
"I X 2 (l-x) ^ 2 (1 + x) 

d<i>i{D,x) _ 2{D~2) 
dx X 

d^^{D,x) _ 2{D - 2) 
dx X 



dx 1 X 2{l-x) 2 (1 -t-x) j ' ' ' 1 X 2{l-x) 

-y^^'^(^3<f2{D,x)-3<fs{D,x) + MD,x)-3MD,x)^ ' ^^^^^ 

d<^2iD, x) _ {D- 2) (^^^^^ _ 2<i,3(^, 3.) + $4(i^, x) - 2$5(i^, x)^ , (172) 

<^i{D,x) - 3^2{D,x) - ^i{D,x) 

<^3{D,x)^<^>^{D,x)^ , (173) 

^2{D,x)+^i{D,x)'j , (174) 

<S>s{D,x)+<i>5{D,x)) , (175) 
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The system is homogeneous: indeed, quite in general the non homogeneous terms are given by 
the Mi's of the subtopologies of the considered graph, obtained by shrinking to a point any of its 
propagator Hnes. When that is done for the parent topology with five propagators of Fig. ([7]), one 
obtains the product of 4 tadpoles; but as the considered graph has two massless propagators, at 
least one massless tadpole is always present in the product; as in the ZJ-dimensional regularization 
massless tadpoles vanish, the product of the 4 tadpoles is always equal to zero - and therefore the 
differential equations are homogeneous. 

By inspection, one sees that $3(1?, x), $5(1?, x) appear in the r.h.s. of Eq.s (|17Hll75p only in the 
combination 

«'3(Z?,x) = $3(^,a:) + $5(^,a;) ; (176) 

the other linearly independent combination of the two Mis, say 

■^^{D,x) = ^:i{D,x)-<i>^{D,x) , (177) 

decouples and can be expressed in terms of the other integrals by means of the trivial 1st order 
differential equation 

d*5(D,x) _ f3(i:>-2) 3(i:>-2)' 
dx ^ ^1 2 (1 - x) ~ 2 (1 + x) 

9 {D-2) 9 {D-2) 
~ j X ^ 2 (1-x) ^ 2 (1 + x) 

As 4'5(D,x) does not enter in the r.h.s. of Eq.s (|17im75p . the 4 linear equations for <I>i(I?,x), 
$2(1?, x), \E'3(Z3,x), and $4(1?, x) can be written as a fourth order equation for $i(x), which will be 
called simply $(!?, x) from now on, and which is therefore equal to 

One obtains for $(_D,x) the following fourth-order ODE 

x3(l - .^)'^^ + x2{l + 5x2 _ „ _ 3x-)Y^^ 

(f<^>{D,x) 



<i>i{D,x)-5'i>2{D,x)-'i>i{D,x) 

'^3{D,x) (178) 



-x|l2 + 6x2 + (£) _ 4)(i3 + 32a;2) + {D - 4)2(1 + 26x2 
+ |l2 - 18x2 + {D- 4)(25 - 2x2) + 8{D - 4)2(2 + 5x2)+ 

+3(15-4)3(1 + 8x2) 



dx^ 



1 d^D,x) 



dx 

4x<{ + 12 + 29{D - 4) + 23{D - 4)2 + 6{D - 4)^ \<^{D, x) = 



(180) 
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9.2. Behaviour of ^{D,x) in the limit x — > 

By inspection, one finds that the most general solution of Eq. (ll80[) can be expanded for x in 
the form 

4 / oo \ 

HD,x)=J2x"'iY^ , (181) 

i=l \n=0 / 

where the values of the 4 exponent at are 



"1=0, 

a2 - (i? - 2) , 

as = -iD-2) , 

a4 = (3I?-7); (182) 



the A};>{D) are the 4 integration constants, and all the other coefficients Ail {D) for n > are 
determined by the differential equation Ea. (|180p . once the integration constants are fixed. 

A qualitative inspection of the integrals which one tries to evaluate by means of the differential 
equation Ea. (|179p shows that it is finite (just finite, not analytic!) for x ^ 0+ and (£> — 2) > 0; that 
is sufficient to rule out from their expression as solutions of the differential equation the terms with 
the behaviour of the third and the fourth exponent (which is negative when D is just above 2). 

In the current the equation for $(!?, x) is homogeneous, the only information one gets out 

is that A'g\D) and A'^\d) are both equal to zero, due to the finiteness for x 0+; by substituting 
the ansatz in Eq. lfTSO]) and dropping A'-q^ {D),A''q^ (D), one finds for <^{D, x) Eq. (fT79l) the 

X expansion 



+ A'fin) (i + + 0(.«) I . (183) 



The expansion depends on the two as yet unspecified integrations constants Aq (D), Aq {D). To fix 
them, one has to provide some independent information, such as the value of the required Feynman 
integral and of its first derivative at a; = 0. Those value can be provided by an explicit conventional 
calculation, say in parameter space, which is in any case much easier than a calculation for non-zero 
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Fig. 8. Four-loop watermelon diagram: a wavy line stands for massless propagator; a solid line, for propagator of 
mass A/. 

values of the variable x 0. That is done explicitly in [133^ and the results are 

,(1) _ 3i^-ll 



4^ 
^0 



8(1? - 2){D -3){D- A)3{2D - 5)(2D - 7)(3L» - 8)(3L» - 10) 

r(i -{D- 4))r(i - 2(D - 4))r2 (i + ^{d ~ 4)) (i - ^{d ~ 4)) 

^ r4 (1 - - 4)) r(l - 3{D - 4)) 

^(2) ^ 2 {2D - 7) 

" 3{D ~ 2)2(1? - 3){D ~ 4)4(37:> - 8){3D - 10) 

r (1 + ^{D - 4)) r (1 - |(z? - 4)) r2(i - (i? ~ 4)) 

^ r2(l-i(D-4))r(l-2(D-4)) 
where the term A^'^ is the value of the vacuum graph in FiglHl 



(184) 



9.3. Expansion around D — 4 and homogeneous equation 

The Laurent's expansion in (D - 4) of ^{D,x) Eq. ()179|) is 

oc 

$(L»,x) = J2 (£'-4)"$(")(a;) , (185) 

as it is known on general grounds that it develops at most a fourth order pole in {D — 4). By 
substituting in Eq. (jl80p one obtains a system of inhomogeneous, chained equations for the coefficients 
$(")(a;) of the expansion in {D — 4); the generic equation reads 



x»$(«)(x) = iv:(")(x) 



(186) 



with 



V = 



6x{2^ 



Sx^) — +48a; 
ax 



(187) 



''In the present case the knowledge of the regularity of the solution at x = 1 does not provide any additional 
information. 
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and 



+ {92x + (16 + 40x^) A _ + 26.^) ^]^^-'^ (x) 



+ (ll6x+ (25 - 2x2) A _ A 3^ ^ 32x3)-^ - (3x^ - g^^) -^U("-i)(a;) 

(188) 

which shows that the equation at a given order n for <i>(")(a;) involves in the inhomogeneous term 
the coefficients ^^'^^x) (and their derivatives) with k < n (obviously ^^''\x) = when k < —4). 
Such a structure calls for an algorithm of solution bottom-up, i.e. starting from the lowest value of 
n (which is n = — 4) and proceeding recursively to the next n + 1 value up to the required order. 
The Eq.s (|186p have all the same associated homogeneous equation, independent of n, 



- ^^)-TT + a;2(l + 5x2)-^ - 6x(2 + x^)-^ + 6(2 - 3x^)4- + 48x 
dx* dx-^ dx'^ dx 



(x) = ; (189) 



once the solutions of Eq. ()189p are known, all the Eq.s (ll86p can be solved by the method of the 
variation of the constants of Euler. 

To our (pleasant) surprise, the solutions of Eq. (|189|l are almost elementary. By trial and error, a 
first solution is found to be 

(j)i{x)=x'^ . (190) 

We then substitute 4'{x) = (^i(x)^(x) in Eq. p89|) . obtaining the following 3rd order equation for the 
derivative of ^ (x) 



(73 ^2 ^ 

x3(l - x^)— r + 3x^(3 -x^)—r+ 6x(l + 2x2)— - 6(5 + 2x^) 
dx-^ dx'^ dx 



1 



C'(x) - , (191) 
(192) 



and find that it admits the solution 

e2{x)^-{l-x'+x'') . 
X'^ 

Substituting ^'(x) = C2{x)x{x) ™ Eq. (|19ip we obtain the following 2nd order equation for x'(x) 

- x2 + x^)-^ + 6x^(2 - x^)4- - 6(2 - 2x4 + x^)] x'{x) - , (193) 
ax^ dx J 

which admits as solution 

' ^-5-2x2 + 5x4 ^^g^^ 



Xsix) = ^(1-2; ) 2 , 4^2 ■ 

x-^ (1 — x^ + x*)^ 

Finally, substituting x'i^) = x'3{x)t{x) in Eq. (|193p . we obtain the equation 



x(l - x'fil -x'+ x*)(5 - 2x2 ^ 5^4) 



dx 



-2(1 - x2)4(15 - 12x2 ^ ^^^4 _,_ 3Q^6 _ 24a,8 20x^ 



t'(x) = , 



(195) 
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which has the solution 



(l-a;2)5 5-2a;2 + 5a;4 



(196) 



By repeated quadratures in x and multipUcations by the previous solutions we obtain the explicit 
analytic expressions of the 4 solutions of Eq. (|189p : the nasty denominators (1 — x^+a;"*) and (5 — 2x2 + 
hx^) disappear in the final results, while the repeated integrations of the terms with denominators 
a:, (1 + x) and (1 — x) generate, almost by definition, HPL's of argument x and weight up to 3. 



(t>2{x) 

03 (a;) 

04 (a;) 



]^{l-x^)-H{0,x)x^ 



(5 + 18x2 + Ux^ + 5a;' 
8a?2 



.8^ 1 

+ -(12 + a;2 - 12a;'')i?(0; x) + 12 x^H{Q, 0; a;) 



(197) 
(198) 



(l + x2)(15 + 182^ 2 + 15^4) , 3(l-a;2)2(5 + x2)(l + 5x2) 

65536 a; 
9(1 -a;4) r 



8192 



131072 a;2 
9x2 , 



H{-l-x)+H{l-x) 



4096 



i7(0,0,-l;x) + iJ(0,0, l;x) 



i?(0,-l;x) + i?(0,l;x) 
The corresponding Wronskian has the remarkably simple expression 

cj)i{x) 02 (a;) 03 (a;) 04 (a;) 

0Ua;) 0'2(x) 0^,(x) 0Ua;) 

W{x) = 

0'/(x) 0i'(x) 0J((x) cj,'i{x) 

c^'i'ix) 0-(x) 0-(x) 0:i"(x) 
in agreement (of course) with the coefficients of the 4*'* and 3'"'^ x-derivative of 0(x) in Eq. (|189p 



(199) 



(1-X2) 



2\3 



(200) 



9.4. Solutions of the chained equations 

With the results established in the previous Section one can use Euler's method of the variation of 
the constants for solving Eq.s ()186|) recursively in n, starting from n = —4. We write Euler's formula 
as 



4 

<i>(")(x) = V 0,(x) 



1 



(n) 



dx' 



W{x' 



-M,{x')K^^Hx') 



(201) 



where the 0i(x) are the solutions of the homogeneous equation given in Ea.s (|190ll99p . the are 
the as yet undetermined integration constants, the Wronskian W{x) can be read from Ea. (|200p . 
the Mi{x) are the minors of the 0j (x) in the determinant Ea. (|200[) . and the K'-^-^x) are the 
inhomogeneous terms of Eq. (|188p . The constants are then fixed by comparing the expansion in 
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a; for a; ^ of Eq. (|20ip with the expansion in {D - 4) for Z? — > 4 of Eq. (|183p . Explicitly, we find, 
for the first two coefficients of the Laurent expansion in Eq. (|185p . 

= -j^^^ , (202) 

(204) 

The full results become quickly too lengthy to be reported explicitly here. 



9.5. Value at x — 1 from Difference Equations 

In this section we will calculate ^{D,x — 1) by using the Difference Equation technique described 
by Laporta in which is a formidable way to get numerical results with high accuracy. We will 
see in this example, the link between Differential Equations and Difference Equations for Feynman 
integral. 

Having set Di = {kf+l),B2 = ((fc2-fci)2 + l), B3 = {{h-k2)^+l),Bi = ih-ksy,!}^ = (p-fc4)^ 
we define 



-2D 



d^ki d^k2 d^ks d^ki 



(205) 



SO that Ic,{D^l) is equal to $(£>, 2:) of Ea. (|f 79p at a; = f up to a known multiplicative factor 

h{D,l) = [4r(f + e)]'^$(L> = 4-2e,a;= f) . (206) 



By combining identities obtained by integration by parts one finds that 15(0, n) satisfies the third- 
order difference equation 



32(n- l)(n-2)(n-3)(n-3D + 5) hiD,n) 



-4(n-2)(n-3) 



f + (39 - 50i:»)n + 27D^ + 50-54 



I2n^ - {38D + 24)n^ + {23D^ + 133D - 84)n 



- lAlD^ + f 34i:> - 24 



h{D,n-2) 



+ {n-D- l){n - 2D + l){2n - 3D){2n - 5D + 4) h{D,n-?,) = 
We will solve this difference equation by using the Laplace ansatz 



(207) 



h{D,n) 



f'v5{D,t) dt , 



(208) 
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giving for vc,[D,t) the fourth-order differential equation 



24(L> + 1)^2 + (18 - 26D)t -ID + 12 



{576{D - l)t^ + (-1081)2 - 420i:i + 648)t^ 



-{A6(f + 38D - 144)< + 71D^ - 284D + 288 



dt 



+t 



{576D - 960)t^ + (-2161)2 + 360D)t^ 



{-18D^ + 190D^ - 496L» + 384)< + 7713^ - 533D^ + 1236D - 960) 



Tt^^iD,t) 



+ {D - 3)(2D - b){2,D - 8){bD - 12) v^{D, t) ^ 



(209) 



We win fook; for the solution of Eq. (|209p in the form of a power series expansions, which, inserted 
in Eq. (|208p and integrated term by term, will provide very accurate values of Ic^{D,n). As the 
convergence is faster for larger n, we will consider large enough values of the index n (see below); 
the repeated use top-down of Eq. p07p {i.e. using it for expressing Ir,[D,n — 3) in terms of the 
75(fc) with A; = n, n — 1, rt — 2) will give the values corresponding to smaller indices, till ^{D, 1) is 
eventually obtained. To go on with this program, initial conditions for v^{D,t) are needed. 

From the definition Eq. (|205p . and introducing spherical coordinates in D-dimension for the loop 
momentum fci, d^ki — kP^^dkidfl{D, ki) one has 



h{D,n) = 



(210) 



iD,kl)^l'-^^UiD,l,iP^h 



(211) 



where r2(_D) is the D-dimensional solid angle, and Ii{D, 1, {jj — fci)^) is the 3-loop (off mass-shell) 
sunrise integral 

/4(A«,(p-fci)2) = 7r-^^/2 / ^""^^ f ^""^^ . 
By the change of variable l/(fc2 + 1) = t, ki — {\ — t)/t, one finds 



(212) 



h{D,n) 



1 



1 - t 



r(f)7o ^ ' -V ^ 

from which one gets the relation between vc,[D, t) and /^{D, {1 — t)/t) 



dt 



dt 



(213) 



(214) 



From that relation we see that we can derive boundary conditions for V5{D,t) in the t ^ 1 limit 
from the expansion of f5{D^k\) in the fci ^ limit, which is easy to obtain. Indeed, only the 
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first denominator of Eq. (l212p depends on fci; expanding for small ki and performing the angular 
integration one gets 



dn{D,ki) 



1 



n{D) (fc2 - + 1 



^2 



dn{D,ki) 
n{D) 

+ kj 



1 



kl - 2fci • fca (fc? - 2ki ■ k2f 
4 



fcf 



{kl + iY i^(fc| + l)3 



+ 



(215) 



(216) 



The above result gives the expansion of f^{D,k\) at kl ~ 0: 

h{D,kl) = /i°^(i^) +/i'^(i^) kl + 0{kt) , 

Note that f5{D, kf) is regular in the origin. 

By inspecting the differential equation p09p one finds that the behaviour at i = 1 of the 4 
independent solution is (1 — t)"S with ai = D/2 — 1, a2 = D/2, as = 0, and — 1; for 
comparison with Eq. (l214p the behaviours = 0, and — 1 are ruled out and the expansion reads 



v,{D,t) = (1 - t)*-i [vf\D) + vi'\D) {l-t)+ 0(1 - tr) 
by comparison with Eq. (|216p {t = 1 corresponds to kf — 0), one obtains 



r(#) 



(217) 



(218) 



The values Ii{D^n) of l4{D,n,p^) at p^ = —1 are therefore required 



h{D,n)=hiD,n,p')=7r 



2^ ^ ^-3D/2 



d"k2 d'^ks d'^ki 



1 



(219) 



■^2 -r J-7 iiJ'3JU'4JU^5 

The problem of evaluating the Ii{D, n) is similar to the original problem of evaluating the Iz{D, n), 
but in fact it is much simpler, as the I^iD, n) involve one less loop and one less propagator. As above, 
by using integration- by-parts identities one finds that Ii{D,n) satisfies the third-order difference 
equation 

6(n - 1) (n- 2) (n-3)/4 (£»,«) 
-{n - 2)(n - 3)(10n -7D~ 10)/4(£', n - 1) 
+ {n - 3)(2n2 + {2D - 18)n - 7D'^ + 29D - 8)/4(L», n ~ 2) 
+ {n - D - l){n - 2D-I- l)(2n - 2,D)Ii{D,n-i) = . 

We solve the difference equation by using again the Laplace ansatz 



h{D,n) 



t"-'^V4{D,t) dt 



(220) 
(221) 



Feynman Diagrams & Differential Equations 49 



where Vi{D^t) satisfies the differential equation 



(fi 

+f(t - l)(36t2 + (6 - 7D)t -9D + 18)— .vJD, t) 



d 



+t{m^ - lADt^ + {-7D^ + 331) - 36)t + 13i:i^ - 611? + T2)—va{D, t) 



dt 



+ {D - 3)(2L> - 5)(3Z) - S)vi{D, t) = 
Following the procedure used above, we write 

1 



h{D,n) = 



h{kl) 



oo 2\Z5/2-l 



dn{D, fca 

n{D) 



h{D,{p-k2 



■h{D,{p-k2f) , 
d^kz d^ki 



Vi{D,t) 



((fc3-fc2)2 + l)(fc4-fc3)2(p-fc4)2 



(222) 

(223) 

(224) 
(225) 
(226) 



At variance with the previous case, the function fi{D, fc|) is not regular for fc2 — 0, as at fc2 = 
the value of the external momentum squared (p — k2f' becomes the threshold of the 2-loop sunrise 
graph associated to I'i{D,p^). But it is not difficult to evaluate analytically l3{D,q^) for generic 
off-shell by using Feynman parameters: 



2T{b-D)T{i-%)T^{%-l 
{D^Af{i-D)T(%) 



D D 

■2F1 (3-i?,2--;-;-q^ 



where 2F1 is the Gauss hypergeometric function. The expansion of I^{D,q^) in q^ = —1 consists of 
the sum of two series. 



h(D,q^)^a^{D) 



l + 0{q^ + l] 



+ bo{D) (q^ + 1) 



2D-5 



1 + 0(92 + 1) 



(227) 



2r(5 - D)T (3 - f ) (f - 1) T{2D - 5) 
(4 - D)^{3 - D)T (fZ? - 3) T{D - 2) ' 
D 



aoiD) - h{D,-l) 



bo{D) = r'l^--ljr{5-2D) . 

Inserting Eq. (|227p into Ea. (|224p . setting q = p — k2 and performing the angular integration over ^2 
by means of the formula (see Eq.(88) of Ref.^S) valid for fc2 ^ 



1 



dn{D,k2) 



-^F(^) 



n{D) J {{p - k2Y + l)N --''^> 2V{N)T {\{D ~ N)) ' 



; 



(228) 
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with N — —(2D — 5), as in the term {q^ 



U{D,kl) = ao{D) 



l + 0{k 



^)2D-5 Eq.JllTl), one obtains 



)r(|-^) 



2r(5-27:>)r(i(3i:'-5)) 



l + 0{h 



(229) 



Using the variable l/(fc2 + l) = t in Eq. (|229p and inserting it in Eq. (|226p one gets the initial condition 
for V4^{D, t) at the singular point t ~ 1 



Vi{D,t) 



r(f) 



T{l-D 



l + 0{l-t) 



2ri5 - 2D)r {^{3D - 5)) 
By inspecting the equation (|222|) one gets that the behaviour at t 



l + 0(l-t) 



ofv4{D,t) is 



(230) 



(231) 



so that for ^ 4 the integral (|22ip is convergent for n > 4. 

All the quantities depending on D are then systematically expanded in (Z) — 4) ; and the series 
are truncated at some fixed number of terms. We solve finally the differential equation (I222p with 
the initial condition (|230|) by a first expansions in series aX t = 1; due to the presence in Eq. (|222p of a 
singular point at t = —1/3, in order to have fast convergence till t — 0, we switch to the subsequent 
series expansions at the intermediate points 1/2, 1/4, 1/8 and 0; then we calculate the integral (|22H) 
for n = 4, 5, 6, 7, 8 by integrating the series term by term By applying repeatedly top-down the 
recurrence relation ([220]) to 14(0,8), h{D,7), Ii{D,&), we obtain h{D,'b) and h{D,4), h{D,3), 
Ii{D, 2) and l4,{D, 1). Those values of /4(_D, n) are used to determine the initial condition for V5{D, t), 
Eq.s (|217|218|216|) . The game must be repeated again for V4{d,t). In fact, we solve the differential 
equation (|209p by expansions in series centered in the points t — 1, 1/2, 1/4, 1/8, 1/16 and (as 
above, this subdivision is due to the presence of a singular point at t — —1/8). By inspecting the 
equation (|209p one gets that the behaviour at t = of V4{d,t) is 



+ (i?)<(-3I3+8)/2 ^ ^(4) p)^(-5I3+12)/2 ^ (232) 



„(4) 



SO that, when D ^ 4, the integral (|208p is surely convergent for n > 5; then we calculate the integral 
P08p for n — 5,6,7,8,9 by integrating the series term by term. By using repeatedly top-down 
the recurrence relation (|220p starting from n = 9, we finally obtain Iz{D, 6), I^{D, 5), I5{D, 4), . . ., 
I^lD,!). By taking into account the normalization (|206p one finds complete agreement with the 



value at a; = 1 of solution of the differential equation computed in the previous section. 



10. Conclusions 

The evaluation of multiloop Feynnian diagrams in the last years has received a strong boost, thanks 
to the ability of turning generic properties of scalar integrals in dimensional regularization into 
tools for computing them. Integration-by-parts, Lorentz invariance, and kinematic symmetries have 
been exploited to establish infinite sets of relations among integrals sharing (partially) common 
integrands. The Laporta algorithm systematizes the by now standard reduction to Master Integrals, 
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that is the algebraic procedure for expressing any Feynman integral as a linear combination of few 
basic integrals with the simplest integrands. 

The completion of the computational task, consisting in the actual evaluation of the Master 
Integrals can be as well afforded by employing the same set of identities. In fact, by combining the 
differentiation of Master Integrals with respect to their Mandelstam invariants, and the reduction of 
the new born integrals, it is possible to derive a system of non- homogeneous first order differential 
equations fulfilled by the Master Integrals themselves. 

Solving such a system of differential equation amounts to the evaluation of the Master Integrals, 
alternatively to their direct parametric integration. 

We have reviewed the method of differential equations by its direct application, trying to follow 
a didactical path. Wc discussed the reduction algorithm plus the general derivation of differential 
equations for Feynman integrals. Successively, the calculation of Master Integrals in the context of 
the evaluation of the one- and two-loop corrections to the photon propagator in QED; whereas, 
in the last two sections, we presented two cases of less trivial differential equations, to show more 
technical aspects related to the solution of homogeneous equations and to the choice of the boundary 
conditions. 

In general, solving a system of first order differential equations for more than one Master Integral 

is equivalent to solving a higher order equation for a single Master Integral. Despite to the lack of 
a theoretical procedure for solving differential equations in the most general case, Euler's variation 
of constants offers a viable procedure. Accordingly, the solution of the non-homogeneous equation 
is obtained by quadrature, using as a kernel the Wronskian of the associated homogeneous equation 
- whose solution can be preliminarily found by constants' variation as well. 

The main achievement is the integration of the differential equation for exact value of the pa- 
rameters (Mandelstam invariants and dimensional-parameter). When that is not possible, one can 
Laurent expand the equation, which then becomes a chained system of equations for the Laurent 
coefficients of the solution, suitable for a bottom-up solving algorithm, starting from the lowest 
Laurent coefficient. 

As a natural feature of Euler's variation of constants, the solution manifests an analytic integral 
representation, generic of transcendental functions: a flexible nested structure of multiple integra- 
tions (or equivalently, iterative summations) which benefits of the shuffle algebras induced by the 
integration-by-parts. Within this framework, the actual efforts required by the computation are 
the finding out of the homogeneous solutions, and the definition of new occurring fiuictions, or- 
dered according to their increasing trascendentality - as required by the iterative fulfilment of the 
non-homogeneous equation. 

Boundary conditions are found by imposing the regularity or the finiteness of the solution at the 
pseudo-thresholds of the Master Integrals. This qualitative information is sufficient for the quan- 
titative determination of the arbitrary integration constants. At the diagrammatic level, boundary 
conditions usually correspond to integrals related to simpler diagrams. 

The use of Differential Eriuation in the; external invariants is a very powerful tool for computing 
Feynman (master) integrals. Dimensional regularization was fundamental for the derivation of the 
differential equations we discussed in this review. 

In principle differential identities for integral functions can be derived whenever it is allowed 
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by the algebra of the integral representation under use - as induced by integration-by-parts. And 
their use is not limited to the perturbative description of Feynman diagrams. Therefore we like to 
conclude by remarking that the use of Differential Equations for integrals' evaluation, is not just a 
technique, but a point of view from which any integral is seen under a new perspective, where there 
appear, explicitly exposed, its analytic structure, its singularities which finally determine its value. 

11. Acknowledgment 

We wish to thank Roberto Bonciani, Thomas Gehrmann, and Ettore Remiddi, for their careful 
reading and comments on the manuscript, but especially for their invaluable collaboration during 
these years. The research of PM was supported by Marie-Curie-EIF under the contract MEIF-CT- 
2006-024178. MA wishes to anknowledge the Institute for Theoretical Physics of Ziirich for the 
kind hospitality while part of this work was performed, with the support of the abovementioned 
Marie-Curie-EIF founds. 



References 

1. A. V. Kotikov, Phys. Lett. B 254 (1991) 158. 

2. E. E. Boos and A. I. Davydychev, Theor. Math. Phys. 89 (1991) 1052 [Teor. Mat. Fiz. 89 (1991) 56]. 

3. A. V. Kotikov, ITF-91-26E 

4. A. V. Kotikov, Phys. Lett. B 259 (1991) 314 

5. F. A. Lunev, Phys. Rev. D 50 (1994) 6589 [arXiv:hep-th/9407174| . 



6. R. Scharf, PhD theis, Wurzburg 1994. 

7. E. Remiddi, Nuovo Cim. A 110 (1997) 1435 [arXiv: hep-th/9711188|. 



8. T. Gehrmann and E. Remiddi, Nucl. Phys. B 580, 485 (2000) arXiv:hep-ph/9912329 

9. T. Gehrmann and E. Remiddi, Nucl. Phys. B 601, 248 (2001) arXiv:hep-ph /0008287 



10. T. Gehrmann and E. Remiddi, Nucl. Phys. B 601, 287 (2001) [arXiv: hep-ph/01 01124 . 

11. C. Ford and D. R. T. Jones, Phys. Lett. B 274 (1992) 409 [Erratum-ibid. B 285 (1992) 399]. 

12. M. Caffo, H. Czyz, S. Laporta and E. Remiddi, Acta Phys. Polon. B 29, 2627 (1998) 
arXiv:hep-th/9807119 . 

13. M. Caffo, H. Czyz, S. Laporta and E. Remiddi, Nuovo Cim. A 111, 365 (1998) a rXiv:hep-th/9805118| . 

14. R. Bonciani, Acta Phys. Po lon. B 30 (1999) 3463. 

15. A. V. Kotikov, arXiv:hep-ph/0102178 

16. M. Czachor and H. Czyz, Acta Phys. Polon. B 32, 3823 (2001) arXiv:hep-ph/0110351). 



17. M. Caffo, H. Czyz and E. Remiddi, Nucl. Instrum. Meth. A 502, 613 (2003) arXiv:hep-ph/0211171 



18. M. Caffo, H. Czyz and E. Remiddi, Nucl. Phys. Proc. Suppl. 116, 422 (2003) arXiv:hep-ph/0211178 



19. P. Mastroha and E. Remiddi, Nucl. Phys. Proc. Suppl. 116, 412 (2003) ar Xiv:hep-ph/0211210 

20. U. Aglietti, arXiv:hep-ph/0408014, 

21. R. Bonciani, Acta Phys. Polon. B 35, 2587 (2004) ,arXiv:hep-ph/0410210| . 

22. V. A. Smirnov, "Feynman integral calculus," Berlin, Germany: Springer (2006) 283 p 

23. M. Caffo, Acta Phys. Polon. B 34, 5357 (2003) 'arXiv:hep-ph/0311052'. 

24. M. Caffo, H. Czyz, A. Grzelinska and E. Remiddi, Nucl. Phys. B 681, 230 (2004) 
arXiv:hep-ph/0312189 . 

25. J. A. Gracey, Phys. Lett. B 277 (1992) 469. 

26. O. V. Tarasov, Phys. Rev. D 54 (1996) 6479 [arXiv:hep -th/96060T8] ; 
Nucl. Phys. Proc. Suppl. 89 (2000) 237 arXiv:hep-ph/0102271 . 

27. O. V. Tarasov, Nucl. Phys. B 502 (1997) 455 arXiv:hep-ph/970 33T9] . 

28. J. Fleischer, A. V. Kotikov and O. L. Veretin, Phys. Lett. B 417 (1998) 163 [arXiv:hep-ph/9707492| . 



Feynman Diagrams & Differential Equations 53 



29. F. A. Berends, A. I. Davydychev and N. I. Ussyukina, Phys. Lett. B 426 (1998) 95 



arXiv:hep-ph/9712209 



30. J. M. Chung and B. K. Chung, Phys. Rev. D 59 (1999) 105014 arXiv:hep-ph/9805432; . 

31. J. Fleischer, A. V. Kotikov and O. L. Veretin, Nucl. Phys. B 547, 343 (1999 ) arXiv:hep-ph/9808242, . 

32. A. I. Davydychev and V. A. Smirnov, Nucl. Phys. B 554 (1999) 391 arXiv :hep^h /9903328| " 

33. J. Fleischer, M. Y. Kalmykov and A. V. Kotikov, arXiv:hep-ph/9905379 

34. J. Fleischer, M. Y. Kalmykov and A. V. Kotikov, Phys. Lett. B 462 (1999) 169 arXiv: hep-ph/9905249| . 

35. J. Fleischer and M. Y. Kalmykov, Comput. Phys. Commun. 128 (2000) 531 arXiv:hep-p h /9907431| 

36. J. Fleischer and M. Y. Kalmykov, Phys. Lett. B 470 (1999) 168 arXiv:hep -ph/9910223^ . 

37. J. Fleischer, O. V. Tarasov and M. Te ntyukov, Nucl. Phys. Proc. Suppl. 89 (2000) 112. 

38. V. A. Smirnov, ■arXiv:hep-ph/0209177 

39. V. A. Smirnov, Nucl. Phys. Proc. Suppl. 116, 417 (2003) [arXiv:hep-ph/0209295 l. 

40. A. Onishchenko and O. Veretin, Phys. Atom. Nucl. 68 (2005) 1405 [Yad. Fiz. 68 (2005) 1461] 
,arXiv:hep-ph/0207091 . 



41. S. P. Martin, Phys. Rev. D 68, 075002 (2003) larXiv:hep-ph/0307101 



42. E. Remiddi, Acta Phys. Polon. B 34 (2003) 5311 arXiv : hep-ph /03 10332 

43. E. Remiddi, Nucl. Phys. Proc. Suppl. 135 (2004) 247. 

44. S. Actis, A. Ferrogha, G. Passarino, M. Passera and S. Uccirati, Nucl. Phys. B 703, 3 (2004) 
'arXiv:hep-ph/0402132 . 

45. T. G. Birthwright, E. W. N. Glover and P. Marquard, JHEP 0409, 042 (2004) [arXiv:hep-ph/0407343| . 

46. B. A. Kniehl and A. V. Kotikov, arXiv:hep-ph/0508238 

47. A. V. Bogdan and R. N. Lee, Nucl. Phys. B 732, 169 (2006) [arXiv:hep-ph/05091811 . 

48. B. A. Kniehl, A. V. K otikov, A. Onishchenko and O. Veretin, Nucl. Phys. B 738, 306 (2006) 
[arXiv:hep-ph/0510235'. 

49. O. V. Tarasov, Phys. Lett. B 638 (2006) 195 |arXiv:hep-ph/0603227| . 

50. C. Anastasiou, T. Gehrmann, C. Oleari, E. Remiddi and J. B. Tausk, Nucl. Phys. B 580, 577 (2000) 
jarXiv:hep-ph/0003261 . 



51. T. Gehrmann and E. Remiddi, Nucl. Phys. Proc. Suppl. 89, 251 (2000) arXiv:hep-ph/0005232 



52. V. A. Smirnov, Phys. Lett. B 500, 330 (2001) arXiv:hep-ph/0011056 . 

53. T. Gehrmann and E. Remiddi, in Proc. of the 5th Inter-national Symposium on Radiative Corrections 
(RADCOR 2000) ed. Howard E. Haber, arXiv:hep-ph/ 0101147 

54. T. Gehrmann, arXiv:h ep-ph/0107090 

55. W. Giele et al., |arXiv: hep-ph/0204316 

56. S. Frixione, Nucl. Phys. Proc. Suppl. 117, 222 (2003) arXiv:hep-ph /0 21 1434| . 

57. J. B. Tausk, Phys. Lett. B 469 (1999) 225 arXiv:hep-ph/9909506 ; 



Z. Bern, L.J. Dixon and A. Ghinculov, Phys. Rev. D 63 (2001) 053007 hep-ph/0010075; ; 

C. Anastasiou, E.W.N. Glover, C. O leari and M.E. T ejeda-Yeomans, Nucl. Phys. B 601 (2001) 

318 [h^-ph/0010212 ; 601 (2001) 347 [hep-ph/00H"094 l; 605 (2001) 486 hep-ph/0101304l; 

E.W.N. Glover, C. Oleari and M.E. Tejeda-Yeomans, Nucl. Phys. 605 (2001) 467 hep-ph/0ld220i; ; 
C. Anastasiou, E.W.N. Glover and M.E. Tejeda-Yeomans, Nucl. Phys. B 629 (2002) 255 



hep-ph/0201274j ; 

E.W.N. Glover and M.E. Tejeda-Yeomans, JHEP 0306 (2003) 033 [hep-ph/0304169| ; 
E.W.N. Glover, JHEP 0404 (2004) 021 hep-ph/04011l9 ; 



Z. Bern, A. De Freitas and L.J. Dixon, JHEP 0109 (2001) 037 hep-ph/0i09078 ; JHEP 0203 (2002) 

018 hep-ph/0201161 ; JHEP 0306 (2003) 028 hep-ph/0304168 ; 

A. De Freitas and Z. Bern, JHEP 0409 (2004) 039 (hep-ph/ 0409007 ; 

Z. Bern, A. De Freitas, L.J. Dixon, A. Ghinculov and H.L. Wong, JHEP 0111 (2001) 031 
hep-ph/0109079 ; 

T. Binoth, E.W.N. Glover, P. Marquard and J.J. van der Bij, JHEP 0205 (2002) 060 [hep-ph/0202266 ; 



54 M. Argeri & P. Mastrolia 



58. L.W. Garland, T. Gehrmann, E.W.N. Glover, A. Koukoutsakis and E. Remiddi, Nucl. Phys. B 627 
(2002) 107 liep-ph/0112081 and 642 (2002) 227 liep-ph/0206067 . T. Gehrmann and E. Remiddi, 
Nucl. Phys. B 640 (2002) 379 jhep-ph/020 7020 . ' 



59. C. Buttar et aL, |arXiv:hep-ph/0604T20l 

28 (1997) 841 arXiv:hep-ph/9610510 

|arXiv:hep-ph/9907494| . 



60. A. I. Davydychev, Acta Phys. Polon. B 28 (1997) 841 arXiv:hep-ph/9610510 



61. C. Anastasiou, E. W. N. Glover and C. Oleari, Nucl. Phys. B 572 (2000) 307 

62. A. V. Kotikov, arXiv:hep-ph/0102177, 



63. G. Passarino, Nucl. Phys. B 619, 257 (2001) arXiv:he p-ph/0108252|; 

G. Passarino and S. Uccirati, Nucl. Phys. B 629, 97 (2002) [arXiv:hep-ph/0112004| ; 



A. Ferroglia, M. Passera, G. Passarino and S. Uccirati, Nucl. Phys. B 680, 199 (2004) 
, arXi vihep-ph /031 1 186 ; 

S. Actis, A. Ferroglia, G. Passarino, M. Passera and S. Uccirati, Nucl. Phys. B 703, 3 (2004) 
farXiv:hep-p h/0402132^ 



G. Passarino and S. Uccirati, Nucl. Phys. B 747, 113 (2006) [arXiv:hep-ph/060312l| 

64. A. V. Kotikov, arifi^hei>ph/0il2347, 

65. S. Laporta, Int. J. Mod. Phys. A 15, 5087 (2000) [arXiv:hep-ph/ 0102033| . 

JarX 



66. S. Laporta, Acta Phys. Polon. B 34 (2003) 5323 [ar Xiv:hep-ph/03 11065 

67. S. Laporta, Phys. Lett. B 549, 115 (2002) [arXiv'he p-ph/0210336' . 

68. A. Ferrogha, M. Passera, G. Passarino and S. Uccirati, Nucl. Phys. B 650 (2003) 162 
arXiv;hep-ph/0209219 . 

69. A. T. Suzuki, E. S. Santos and A. G. M. Schmidt, J. Phys. A 36 (20 03) 11859 la rXiv:hep-ph/0309080| . 

70. V. A. Smirnov, Nucl. Phys. Proc. Suppl. 135, 252 (2004) [arXiv: hep-ph/0406052]^ 



71. T. Binoth and G. Heinrich, Nucl. Phys. B 585 (2000) 741 'arXiv:hep-ph/0004013^ 

72. T. Binoth and G. Heinrich, Nucl. Phys. B 680 (2004) 375 arXiv:hep-ph/0305234 

73. G. Duplancic and B. Nizic, Eur. Phys. J. C 35, 105 (2004) ]ar'Xiv:hep- ph/0303184 

74. M. Czakon, arXiv:hep-ph/0511200 



75. C. Anastasiou and A. Daleo, arXiv:hep-ph/0511176 



76. V. A. Smirnov, Nucl. Phys. Proc. Suppl. 157, 131 (2006) [arXiv:hep-ph/0601268l . 

77. S. Weinzierl, arXiv:hep-ph/0604068 

78. A. G. Grozin, Int. J. Mod. Phys. A 19 (2004) 473 arXi v:hep-ph/0307297| . 

79. V. A. Smirnov, Springer Tracts Mod. Phys. 211 (2004) 1. 

80. R. Bonci ani, P. Mastro ha and E. Remiddi, Nucl. Phys. B 661, 289 (2003) [Erratum-ibid. B 702, 359 
(2004)] arXiv:hep-ph/0301170 . 

81. R. Bonciani, P. Mastroha and E. Remiddi, Nucl. Phys. B 676, 399 (2004) arXiv:hep-ph /0307295| . 

82. R. Bonciani, Nucl. Phys. Proc. Suppl. 152, 168 (2006) arXiv:hep-ph/0410 092|. 

83. R. Bonciani, P. Mastrolia and E. Remiddi, Nucl. Phys. B 690, 138 (2004) arXiv:hep-ph/0311145l. 

84. W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, T. Leineweber, P. Mastrolia and E. Remiddi, 
Nucl. Phys. B 706, 245 (2005) arXiv:hep-ph/0406046 . 

85. W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, T. Leineweber, P. Mastrolia and E. Remiddi, 
Nucl. Phys. B 712, 229 (2005) [arXiv :hep-ph/0412259 . 

86. W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, T. Leineweber and E. Remiddi, Nucl. Phys. 
B 723, 91 (2005) arXiv:hep-ph/0504190 . 

87. W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, T. Leineweber, P. Mastrolia and E. Remiddi, 
Phys. Rev. Lett. 95 (2005) 261802 ar Xiv:hep-ph/05093 4i; . 

88. A. V. Kotikov, A. V. Lipatov, G. Parente and N. P. Zotov, Eur. Phys. J. C 26 (2002) 51 



arXiv:hep-ph/0107135 



89. A. V. Kotikov and L. N. Lipatov, Nucl. Phys. B 661 (2003) 19 [Erratum-ibid. B 685 (2004) 405] 
;arXiv:hep-ph/0208220;. 

90. F. Jegerlehner and M. Y. Kalmykov, Nucl. Phys. B 676 (2004) 365 [arXiv:hep-ph/0308216l . 



Feynman Diagrams & Differential Equations 55 



91. C. Anastasiou, L. J. Dixon, K. Melnikov and F. Petriello, Phys. Rev. D 69 (2004) 094008 
[aSiv:hep-ph/0312266 . 

92. K7 G. Chetyrkin, J. H. Kuhn, P. Mastrolia and C. Sturm, Eur. Phys. J. C 40 (2005) 361 
^arXiv:hep-ph/0412055 . 

93. YTSchroder and A. Vuorinen, JHEP 0506 , 051 (2005) [arXiv:hep-ph/050"3209] . 

94. M. Faisst, P. Maierhoefer and C. Sturm, arXiv:hep-ph/0611244' 

95. M. Neubert, Phys. Kept. 245 (1994) 259 arXiv:hep-ph/9306320 . 

96. A. Juste et al, arXiv:hep-ph/0601112 

97. W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, T. Leineweber, P. Mastrolia and E. Remiddi, 
PoS HEP2005, 229 (2006) arXiv:hep-ph/0601207 . 

98. W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, T. Leineweber, P. Mastrolia and E. Remiddi, 
arXiv:hep-ph/0604031 

99. R7 Bonciani, arXiv:hep-ph/0607037 

100. U. Aglietti and R. Bonciani, NucL Phys. B 668, 3 (2003) [arXiv:hep-ph/0304028' . 

101. U. Aglietti and R. Bonciani, Nucl. Phys. B 698, 277 (2004) arXiv:hep-ph/0401193 . 

102. U. Aglietti, R. Bonciani, L. Grassi and E. Remiddi, arXiv:0705.2616 [hep-ph]. 

103. P. N. Maher, L. Durand and K. Riesselmann, Phys. Rev. D 48 (1993) 1061 [Erratum-ibid. D 52 (1995) 
553] [ajXiv?he p-ph/9303233 . 

104. U. Aglietti, R. Bonciani, G. Degrassi and A. Vicini, Phys. Lett. B 595, 432 (2004) 
^iv:hep-ph/0404071 . 

105. U. Aglietti, R. Bonciani, G. Degrassi and A. Vicini, Phys. Lett. B 600, 57 (2004) 
arXiv:hep-ph/0407162^. 

106. W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, P. Mastrolia and E. Remiddi, Phys. Rev. D 
72, 096002 (2005) 'arXiv:hep-ph/0508254 . 

107. C. Anastasiou, S. Beerli, S. Bucherer, A. Daleo and Z. K unszt, |arXiv :hep-ph/0611236[ 

108. U. Aglietti, R. Bonciani, G. Degrassi and A. Vicini, arXiv:hep-ph/061 1266[ 

109. V. A. Smirnov, Phys. Lett. B 524, 129 (2002) arXiv:hep-ph/0111160 ; . 

110. R. Bonciani, A. Ferrogha, P. Mastroha, E. Remiddi and J. J. v an der Bij, Nucl. Phys. B 681, 261 
(2004) [Erratum-ibid. B 702, 364 (2004)] arXiv:hep-ph/0310333; . 

111. R. Bonciani, A. Ferrogha, P. Mastroha, E. Remiddi and J. J. van der Bij, Nucl. Phys. B 701, 121 

(2004) arXiv:hep-ph/0405275 . 

112. G. Heinrich and V. A. Smirnov, Phys. Lett. B 598, 55 (2004) [arXiv:hep-ph /0406053 I . 

113. R. Bonciani, A. Ferro gha, P. Mastroha, E. Remiddi and J. J. van der Bij, Nucl. Phys. B 716, 280 

(2005) arXiv:hep-ph/0411321 . 

114. M. Czakon, J. Gluza and T. Riemann, Phys. Rev. D 71, 073009 (2005) [arXiv:hep-ph/0412164 . 



115. R. Bonciani and A. Ferrogha, Phys. Rev. D 72, 056004 (2005) arXiv:hep-ph /050704 7] . 

116. R. Bonciani and A. Ferrogha, Nucl. Phys. Proc. Suppl. 157, 11 (2006) ;arXiv:he p-ph/0601246" . 

117. M. Czakon, J . Gluza, K. Kajda and T. Riemann, Nucl. Phys. Proc. Suppl. 157, 16 (2006) 
'arXiv:hep-ph/0602102 . 

118. D. Seidel, Phys. Rev. D 70, 094038 (2004) tarXiv: hep-ph/0403185| . 

119. K. Melnikov and A. Mitov, Phys. Lett. B 620, 69 (2005) [arXiv:h ep-ph/0505097 . 

120. H. M. Asatrian, T. Ewerth, A. Ferrogha, P. Gambino and C. Greub, Nucl. Phys. B 762 (2007) 212 



arXiv:hep-ph/0607316] . 



121. F. De Fazio, arXiv:hep-ph/0010007, 

122. F. Jegerlehner, M. Y. Kalmykov and O. Veretin, Nucl. Phys. B 658 (2003) 49 [arXiv:hep-ph/0212319l . 

123. M. Awramik, M. Cza kon, A. Freitas and G. Weiglein, Phys. Rev. Lett. 93, 201805 (2004) 
WXiv:hep-p h/04073T7 ^. 

124. M. Awramik, M. Czakon, A. Freitas and G. Weiglein, Nucl. Phys. Proc. Suppl. 135, 119 (2004) 
arXiv:hep-ph/0408207 . 

125. A. Freitas, M. Awramik and M. Czakon, ,arXiv:hep-ph/050715"9} 



56 M. Argeri & P. Mastrolia 



126. M. Awramik, M. Czakon and A. Freitas, ' arXiv:hep-ph /06053391 

127. D. J. Broadhurst, arXiv:hep-th/9806174 ^ 

128. A. I. Davydychev and R. Delbourgo, Acta Phys. Polon. B 29, 2891 (1998) arXiv:hep"-th/9806248' . 

129. Remiddi, E. and Vermaseren, J. A. M., Int. J. Mod. Phys. A 15 (2000) 725-754 he p-ph/9905 237 . 
Gehrmann, T. and Remiddi, E., Comput. Phys. Co mmun. 141 (2001) 296 |(hep -ph/0107173| ; C omput. 
Phys. Commun. 144 (2002) 200 hep-ph/0111255 ; Nucl. Phys. B 640 (2002) 379~ [hep-ph/0207020| . 

130. D. Maitre, Comput. Phys. Commun. 174 (2006) 222 arXiv:hep-ph/0507152 . 

131. D. Maitre, arXiv:hep-ph/0703052, 



132. Moch, S. and Uwer, P. and Weinzierl, S., J. Math. Phys. 43 (2002) 3363 hep-ph/01 10083 



133. Weinzierl, S., Comput. Phys. Commu n. 145 (2002) 357 ma th-ph/02010111 

134. Blumlein, J., (2003) hep-ph/0311046 

135. M. Y. Kalmykov, Nucl. Phys. Proc. Suppl. 135, 280 (2004) [arXiv:hep-th /0406269] 



136. M. Y. Kalmykov, B. F. L. Ward and S. Yost, arXiv:hep-th/0612240 

137. G. V. Dunne, JHEP 0402, 013 (2004) prXiv: hep-th/0311167 



138. S. P. Martin, Phys. Rev. D 70, 016005 (2004) [arXiv: hep-ph/0312092| . 

139. S. P. Martin and D. G. Robertson, Comput. Phys. Commun. 174, 133 (2006) arXiv: hep-ph/0501 132 

140. G. V. Dunne and M^ Krasnansky, JHEP 0604 (2006) 020 arXiv:hep-th/0602216, . 

141. J. A. Gracey, arXiv:hep -th/0605037 ' 

142. G. Kallen, A. Sabry, Da. Mat. Fys. Medd. 29, No.l7 (1955) 1. 

143. A. Djouadi, P.Gambino Phys. Rev. D49 (1994) 3499 

144. D.J. Broadhurst, J. Fleischer, O.V. Tarasov, Z. fur Physik C60 (1993) 287^ 

145. M. Argeri, P. MastroUa and E. Remiddi, Nucl. Phys. B 631 ( 2002) 388 arXiv:hep-ph/0202123l . 

146. P. Mastroha and E. Remiddi, Nucl. Phys. B 657 (2003) 397 jarXiv:hep-ph/0211451|. 

147. S. Laporta, P. Mastrolia and E. Remiddi, Nucl. Phys. B 688 (2004) 165 'arXiv:hep- ph/0311255| . 

148. S. Laporta and E. Remiddi, Nucl. Phys. B 704 (2005) 349 arXiv:hep-ph/0 406160 . 

149. Chetyrkin, K. G. and Tkachov, F. V., Nucl. Phys. B 192 (1981) 159. 
Tkachov, F. V., Phys. Lett. B 100 (1981) 65. 

150. S. Laporta and E. Remiddi, Phys. Lett. B 379 (1996) 283 [arXiv:hep-ph/9602417l . 

151. P. MastroUa and E. Remiddi, Nucl. Phys. Proc. Suppl. 89 (2000) 76. 

152. SOLVE, by E. Remi ddi, available from the author. 
SYS, by S. Laporta,!^ . 

IdSolver, by M. Czakon. 

AIR, by C. Anastasiou and A. Lazopoulos, JHEP 407 (2004) 04 6 [arXiv:hep-ph/0404258| . 

153. P. A. Baikov, Phys. Lett. B 634 (2006) 325 arXivThep-ph/0507053 . 

154. A. V. Smirnov and V. A. Smirnov, JHEP 0601 (2006) 001 arXiv:hep-lat /0509187| . 

155. R. J. Eden, P. V. LandshofF, D. I. Olive, and J. C. Polkinghorne "The Analytic S-Matrix," Cambridge 
University Press (1966) 287 p 

156. t Hooft, G. and Veltman, M. J. G., Nucl. Phys.", B 44 (1972) 189-213. 

BoUini, C. G. and Giambiagi, J. J., Phys. Lett. B 40 (1972) 566-568; Nuovo Cim. B 12 (1972) 20-25. 

Ashmore, J. F., Lett. Nuovo Cim. 4 (1972) 289-290. 

Cicuta, G. M. and Montaldi, E., Nuovo Cim. Lett. 4 (1972) 329-332. 

Gastmans, R. and Meuldermans, R., Nucl. Phys. B 63 (1973) 277-284. 

Collins, J. C, Renormalization, "Cambridge University Press", (1987). 



